The choroid, a vascular structure between the retina and sclera, is primarily responsible for oxygenation and metabolic activity of the retinal pigment epithelium layer (RPE) and the outer retina. Its main substructures are (i) the Bruch’s membrane; (ii) the choriocapillaris, a capillary network adjacent to the Bruch’s membrane; (iii) the choroidal stroma, layers with intermediate and large vessels that drain the choriocapillaris; (iv) the suprachoroid, a transitional zone between the choroid and sclera ( thick) composed of numerous compactly arranged lamellae of melanocytes, fibroblasts, collagen bundles, and nerve fibers (for more details, see Refs. 12.–3). Long-term variations of choroidal thickness occur in the process of infantile ocular development and healthy aging of adults4,5 as well as in connection with retinal diseases, which can affect the choroid.6,7 In the literature, short-term (even diurnal) thickness variations have been reported as well, cf. Ref. 2, pp. 154 ff.; Ref. 8. Differences in choroidal thickness have also been observed in relation to different levels of myopia.9 Consequently, there is a considerable interest in accurate measurements of choroidal thickness.
Historically, the standard procedure for choroid imaging was indocyanine green angiography, which allows for two-dimensional (2-D) imaging of the choroid pattern and detection of leakage and vessel wall abnormalities.10 Although optical coherence tomography (OCT) has become a standard tool for noninvasive imaging of the retinal layers since its introduction in 1991, its application for the investigation of the choroid has largely been unsuccessful for more than a decade. In OCT images, the (visual as well as automated) identification of the highly reflective RPE/choroid boundary is easy, but the detection of the outer choroid boundary (OCB) is a much more challenging task for several reasons. First, the position of the choroid posterior to the highly scattering layers of the cellular/subcellular structures of the outer retina causes a relatively weak signal from the region as a whole. Second, the observed reflectivity jumps are weak and strongly affected by noise. Only recent advances in OCT imaging, namely real-time eye-tracking (allowing for averaging of large numbers of A-scans), enhanced depth imaging (EDI SD-OCT) (as introduced in Ref. 11), swept-source OCT (SS-OCT)12 or additional use of polarization information13 led to high-quality in vivo visualizations of the choroidal vasculature, providing reliable depth information for choroidal thickness measurements. A third, principal difficulty results from the heterogeneous structure of the suprachoroidal layer, which consists of tightly interlaced components of both low and high reflectivity (melanocytes resp. collagen bands and fibroblasts).3 As a consequence, it may be questioned whether the border between the suprachoroid and the collagen bundles of the sclera, which is histologically well-defined, is always clearly resolved by OCT imaging. Instead, the positive contrast jump visible in the data, which seems to be related to regions with sufficiently high concentrations of melanocytes, is used for the operational definition of the OCB in the present study (as well as in the vast majority of the literature), thus providing a reliable inner approximation of the histological choroid/sclera border.
The built-in segmentation routines in OCT devices are not suitable for scientific measurements as they are designed mainly for diagnostical purposes and, in general, are poorly documented (if at all). For the automated segmentation of retinal layers in OCT images, in particular for the detection of the RPE/choroid boundary, a large variety of scientifically documented methods is available by now. These include denoising by diffusion and subsequent edge detection,14 denoising and simultaneous edge detection by variational methods,15 active contour methods,16 graph-theoretical approaches,1718.19.–20 and classification by support vector machines.21 Furthermore, in Ref. 22, a statistical segmentation method that is essentially based on the application of a neural network was proposed. Note that in almost all aforementioned cases a preceding denoising step is included even when not explicitly mentioned.
In contrast to this situation, the majority of studies concerned with the determination of choroidal thickness rely on manual or semi-manual detection of the OCB.45.6.–7,9,12,23188.8.131.52.184.108.40.206.32.33.34.–35 Unfortunately, manual delineation procedures are quite time-consuming and show considerable intra- and interobserver variability even when based on EDI SD-OCT or SS-OCT data.28,36,37 Therefore, an automated detection of the OCB is highly desirable but has been realized only in a small number of studies until now, applying successful approaches for retinal layer segmentation to OCB detection. The studies36,3839.40.41.–42 pursue graph-theoretical approaches while43 adapted their statistical segmentation technique from their previous work.22
In the present study, we establish a novel method for the automated detection of the OCB within SD-OCT image data. In difference to the studies mentioned above, our approach is based on an image model within the space of functions of bounded variation (BV) and combines the application of quadratic measure filters44220.127.116.11.–49 with elementary feature recognition. To the best of the authors’ knowledge, the method has not been used in the context of medical image processing as yet. It applies to the detection of retinal layer boundaries as well and is particularly suitable for OCT data generated without special imaging modes (as EDI or SS-OCT) and moderate line averaging. To demonstrate the capability of the new approach, only such data (acquired without EDI module, 10 times averaged) were used in the present study. Note that a considerable number of OCT studies employ solely this low-grade image quality (certainly, the method presented applies to data of higher quality just as well, as will be shown in a subsequent publication). Within the B-scans, three layer boundaries were automatically recognized: the inner limiting membrane (ILM) boundary, the RPE/choroid boundary and the OCB. Based on these segmentations, we provide an automated determination of the central fovea region.50,51 Automated measurements of mean choroidal thickness are given for this and two adjacent 1 mm regions, as defined in more detail in Sec. 2.6. The quality of the automated procedures will be assessed by comparison with manual segmentations performed by five trained graders and subsequent thickness derivations. Moreover, we shortly discuss whether the obtained OCB segmentations correspond to the histological structure of the transitional zone between the choroid and sclera, cf. Ref. 2, p. 147; Ref. 3. Our study is based on data from 50 children aged 8 to 13 years, which have been collected in the framework of the LIFE Child study at the Leipzig University.
The outline of the paper is as follows. In Sec. 2, we describe material and methods, particularly providing a detailed description of the automated segmentation method. Section 3 is devoted to the presentation and discussion of the results, Sec. 4 presents our conclusions. For convenience of the reader, the mathematical background of the presented method is briefly summarized in the Appendix. Detailed information about the OCT data is contained in Table 1.
OCT data. Right eye imaged; size of B-scans is 496×512 px; pixel depth is 3.87 μm.
|Subject no.||Fov. center at B-scan (no.)||Fovea pos. (px)||Gender||Age (y)||Mean ax. length (mm)||Scan focus (dpt)||Mean corneal rad. (mm)||Pixel width (μm)|
Material and Methods
Selection of Probands within the Framework of the LIFE Child Study
The LIFE Child study is a population-based longitudinal cohort study conducted in Leipzig, a city in central Germany with more than 500,000 inhabitants. With recruitment age ranging between the 24th week of gestation and 16 years of age and annual follow-ups, the study combines a cross-sectional and a longitudinal design. Participants stem from the city of Leipzig and its close proximity. Recruitment started in 2011 and is planned to continue until 2021. By the end of 2015, more than 3000 children have already participated. Study participants are recruited via advertisement at different institutions such as university hospitals, local clinics, public health centers, kindergartens, schools, and partner study centers. The study program covers the collection of biological samples, examination of body functions and assessment of skills, traits, and habits (for more details, see Ref. 52). For the present work, OCT data from 50 children (25 male, 25 female) aged 8 to 13 have been randomly selected. The use of data was approved by the institutional review board of the Leipzig Research Center for Civilization Diseases (LIFE).
The LIFE Child study follows the tenets of the Declaration of Helsinki. Its outline was approved by the Ethics Committee of the Leipzig University (Reg. No. 264-10-19042010). The study has been registered with the trial number NCT02550236. The parents of the participating children were informed about the study program, the long-term use of data, potential risks of participation and the right to withdraw from the study. Written informed consent is obtained from the parents of all individual participants included in the study.
Acquisition of SD-OCT Image Data
Chorio-retinal images were obtained with a commercially available spectral-domain (SD) OCT (Spectralis HRA + OCT, equipped with camera head with serial number 04514, software modules Heidelberg Eye Explorer 18.104.22.168, Aquisition Module 22.214.171.124 and Viewing Module 126.96.36.199; Heidelberg Engineering, Heidelberg, Germany). The device uses a super luminescent diode with a central wavelength of 870 nm and acquires 40,000 A-scans/s. The scan depth in tissue is 1.9 mm with 496 pixels per A-scan, resulting in a pixel depth of while the axial optical resolution amounts to , cf. Ref. 53, p. 273. For each child, volume scans of the macular area with a field size of were acquired. Each volume scan contains 97 equally spaced averaged B-scans consisting of . Additionally, the device utilizes a confocal scanning laser ophthalmoscope at a wavelength of 815 nm to capture an overview image of the retina. Herein points are mapped to track eye movements using the overview image as a reference. This built-in function for real-time eye-tracking was active in order to obtain an average of 10 B-scans per line. EDI module was not used for the current analysis. For all children, both eyes were measured but the subsequent analysis was carried out on right eye data only. Corneal radii were measured with the Lenstar instrument (LS900 Biometer, Haag-Streit, Wedel, Germany). The median of three to six corneal biometry measurements was used for calculating the mean anterior corneal radius (included in Table 1) by averaging the steep and flat anterior meridians. Axial length was obtained as the mean of three to six measurements per eye, performed by the above mentioned Lenstar instrument. Measurement precision for corneal radii and axial length is about 0.01 mm, cf. Ref. 54, p. 13.
For the purpose of segmentation, exclusively the raw OCT data were exported and used.55 Additionally, for convenience, a horizontal strip of with zero values was joined with each original B-scan. Using the results of the ILM and RPE/choroid boundary detection, from every volume scan, a B-scan defining the foveal center has been singled out (see Sec. 2.5). These scans have been further used for choroidal thickness measurements (see Table 1). A typical B-scan (Table 1, line 45) is shown in Fig. 1(a).
Description and Implementation of the Edge Detection Method
Basic approach. In order to detect a layer boundary within the OCT data , the following five steps were performed: (1) smoothing of by calculation of , (2) calculation of discretized edge detectors , (3) evaluation of the edge detectors by columnwise local maxima search (testing of the alternative expressed in Eq. (34), see Appendix), (4) preliminary estimation of the layer position according to its definition, and (5) final designation of the position by calculation of barycenters and cubic spline interpolation. In the following, the different steps will be reported in full detail.
Step 1. Convolution. The convolution with the scaled Gauss kernel is calculated as the numerical solution of the heat equation on the domain with initial values , cf. Ref. 57, p. 47, Theorem 1. Specifying , the edge length of a pixel is . Note that, at the same time, is the density of the 2-D normal distribution with mean value and standard deviation . Consequently, the region of coincides with the ball centered in the origin and radius for , where .
Step 2. Discretized edge detectors. Assuming that is subdivided into pixels , , , the test functions , are specified throughAppendix A.3, six different edge detectors , , , , , and (“quadratic measure filters") are defined as discretizations of the integrals
Steps 3–5. Recognition of the ILM boundary. The vertical column of pixels (A-scan) at the position is denoted by . For each of the three detectors , within the sets , all local maxima are determined. Their positions are stored within sets . Next, every marked pixel within is replaced by the vertical three-pixel feature . Then the ILM boundary will be preliminary defined as the topmost adjacent feature within of sufficiently large size, marking a positive contrast jump. Denoting this feature by , for the three edge detectors, the barycenters over are calculated, and the mean of the three quantities is formed. The final position is obtained by cubic spline interpolation of the refined data.
Steps 3–5. Recognition of the RPE/choroid boundary. For each of the six detectors , the local maxima over the sets are determined. Their positions are stored within sets . Every marked pixel within is replaced by a cross-shaped feature . Then the RPE/choroid boundary will be preliminary defined as the lowermost adjacent feature within of sufficiently large size, marking a negative contrast jump at the same time. Denoting this feature by , for every edge detector, the barycenters over are calculated. Forming the mean of the six quantities, a refined boundary position is obtained. Again, the final position is found by cubic spline interpolation of these data.
Steps 3–5. Detection of the OCB. Within the usual visualization of the B-scans, this boundary can be visually recognized as a fairly weak positive contrast jump below the RPE/choroid boundary [see Fig. 1(a)]. One may observe that this jump corresponds to a (possibly not connected) feature within the normalized gradient field of where the second component of the gradient vector is nonnegative [see Figs. 1(c) and 1(d)]. Consequently, the maximization of the edge detector will be restricted to a set , which is defined as the union of topmost connected features below the RPE/choroid boundary withFig. 2(a)]. This position is refined in two steps: after calculating the barycenters of over , the final position is obtained again by cubic spline interpolation.
Remarks about the implementation. The complete edge detection procedure was implemented as a series of MATLAB® tools. It was performed with MATLAB® 188.8.131.52421 (R2014b) and required the Image Processing and Signal Processing toolboxes, as documented in Refs. 5859.–60. For the numerical solution of the heat equation, after the introduction of a spatial mesh whose nodes are the centers of the pixels , the ADI method proposed by Peaceman/Rachford was employed, cf. Ref. 61, p. 412 f. Note that the missing smoothness of the initial values does not influence the calculations. For the calculation of and feature recognition, the imfilter and bwlabel procedures were used while the maxima search was realized with the aid of the findpeaks procedure. The typical CPU time for the segmentation of a single B-scan amounts to 7.5 s (steps 2 to 5). No particular attempts for tuning have been made.
Validation of the Computer-Generated Results
The accuracy of the automated detection of the boundary layer positions has been validated by comparison with manual segmentation. To this end, within each of the 50 B-scans from Table 1, the ILM boundary, the RPE/choroid boundary, and the OCB were manually delineated by five trained graders (Anna Xenia Bestehorn, Mike Francke, Anja Schirmer, Sophia Scheibe and Beatrice Zimmerling) independently from each other, employing the usual visualization of the OCT data , cf. Ref. 55, p. 11. The graders were allowed to enhance the contrast within the images if necessary. The manually obtained boundary positions will be denoted by , , and , respectively. In Sec. 3.1, the comparison between manually and automatically obtained data is described in detail.
Definition of the Foveal Region within the B-Scans
Within each volume, the segmentation procedure was applied to the B-scans # 42 … # 57. Subsequently, from each volume the minimal distance of the ILM and RPE/choroid boundaries, compared over the A-scans with numbers , was extracted together with its position, thus defining the foveal center. The respective scan was singled out for further analysis (see Table 1).
In order to define within each of these scans the central fovea region, we rely on a fovea model introduced in Ref. 50. Following the procedure described there, a model function for the fovea shape was generated from the automatically extracted ILM and RPE boundary data and the position of the foveal center. The extent of the fovea bowl left and right from the center is given by two radii and that can be determined from the model (see Fig. 3). As pointed out in Ref. 51, p. 5, instead of using of the highest point of the foveal rim for this purpose, it is more advantageous to solve the equation
Three choroid measurement regions were defined: (1) the central foveal zone, which is the region between and , (2) the 1-mm region left of the central foveal zone (corresponding to temporal retina for the right eyes examined), and (3) the 1-mm region right of the foveal zone (referring to the nasal retina). As discussed in Ref. 51 in detail, one may observe a strong variability of foveal geometry between individuals, requiring the determination of a measurement region of variable width [thus the use of the standard 1 mm circle, which has been introduced in the context of fundus photography and is widely used in retinal thickness measurements (see Ref. 62, Chapter 18; Ref. 63), would be completely inappropriate]. Outside of the foveal rim, the individual shape variation is far less prominent, allowing for the definition of measurement regions of fixed dimensions. In order to obtain the magnification factor that is directly used to calculate the pixel width in mm, we employed the formula published in Ref. 64, p. 649,
Definition and Measurement of Choroidal Thickness
The (pointwise) choroidal thickness is understood as the vertical distance from the hyperreflective line of the RPE/Bruch’s membrane boundary to the positive contrast jump below this curve, which operationally defines the OCB as mentioned above. Within the three foveal regions defined in Sec. 2.5, the mean choroidal thicknesses are defined as
Histological Examination of the Transitional Zone between Choroid and Sclera
After enucleation, human donor eyes were fixed in 4% paraformaldehyde phosphate-buffered solution (PBS) for more than 24 h. Choroidal and scleral tissues were stored in PBS and prepared for light microscopy. Layers of choroidal stroma (Sattler’s and Haller’s layer) were partly removed to get a plane view onto the inner surface of the outer choroid. Plane view and cross-sectional images were obtained using a digital reflected-light microscope (VHX-600, Keyence Deutschland GmbH, Neu-Isenburg, Germany). Use of human tissue was approved by the Ethics Committee of the Leipzig University.
Results and Discussion
Results of the Automated ILM, RPE, and OCB Boundary Detection
In Table 2, we display the results of automated ILM and RPE detection. For every B-scan examined, the pointwise mean and the pointwise standard deviation of the five manual delineations is calculated as , , , and , respectively. With the aid of these quantities, for every B-scan, the averaged standard deviations
Results of automated ILM and RPE detection. Entries defined in Eqs. (11)–(16).
|No.||ILM‾σ (μm)||ILM¯-|Δ|am (μm)||ILM¯-Δam (μm)||RPE‾σ (μm)||RPE¯-|Δ|am (μm)||RPE¯-Δam (μm)|
The results for the OCB are presented in Table 3. Again, for every B-scan examined, the pointwise mean and the pointwise standard deviation of the five manual delineations is calculated as and . Then for every B-scan the averaged standard deviation
Results of automated OCB detection. Entries defined in Eqs. (17)–(20).
|No.||OCB‾σ (μm)||OCB¯-|Δ|am (μm)||OCB¯-Δam (μm)||Perc. (%)||OCB¯-|Δ|a1 (μm)||OCB¯-|Δ|a2 (μm)||OCB¯-|Δ|a3 (μm)||OCB¯-|Δ|a4 (μm)||OCB¯-|Δ|a5|
Results of the Automated Choroidal Thickness Measurements
In Table 4, the choroidal thickness measurements for the central foveal region are presented. For comparison, we determined for every B-scan from the five manual delineations and the thicknesses2.5. For every B-scan, the five values obtained from the manual delineations are tabulated in columns 1 to 5. Columns 6 and 7 contain the mean Table 5 contains the choroidal thickness measurements for the left and right 1-mm regions. Analogous to the computations before, we calculated for every B-scan the thicknesses (left 1-mm region) and (right 1-mm region), , and tabulate for every B-scan the mean
Choroid thickness, central foveal region. Entries defined in Eqs. (9), (21)–(24). Automated measurements in boldface.
|No.||CHC1 (μm)||CHC2 (μm)||CHC3 (μm)||CHC4 (μm)||CHC5 (μm)||CHCm (μm)||CHCσ (μm)||CHCa (μm)||CHC-Δam (μm)|
Choroid thickness, left and right outer 1 mm regions. Entries defined in Eqs. (8), (10), (25)–(30). Automated measurements in boldface.
|No.||CHLm (μm)||CHLσ (μm)||CHLa (μm)||CHL-Δam (μm)||CHRm (μm)||CHRσ (μm)||CHRa (μm)||CHR-Δam (μm)|
Discussion of the Results
General comparison with other methods from the literature. When comparing different methods for OCB segmentation and choroidal thickness measurement, a basic methodical problem is caused by the large variation of quality of the underlying image data [Ref. 36, p. 2869]. The majority of the existing studies is based on data generated by special imaging modes and/or high line averaging (see Table 6). Another problem addressed in the literature is the repeatability of the measurements, which extends to the applied segmentation methods.
Properties of OCT data used in the cited references.
|Reference||Imag. method||Device||Averaging||(Semi-)manual OCB segm.||Automated OCB segm.|
|Manjunath et al.4||SD-OCT||Carl Zeiss Meditec, Cirrus HD-OCT||20||•|
|Manjunath et al.7||SD-OCT||Carl Zeiss Meditec, Cirrus HD-OCT||20||•|
|Hu et al.39||SD-OCT||Heidelberg, Spectralis||9||•|
|Zhang et al.41||SD-OCT||Carl Zeiss Meditec, Cirrus HD-OCT||none||•|
|Zhang et al.42||SD-OCT||Topcon, 3D-OCT 1000 / Topcon, 3D-OCT 2000||none||•|
|Present study||SD-OCT||Heidelberg, Spectralis||10||•|
|Cho et al.23||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Dhoot et al.6||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Ding et al.24||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Flores-Moreno et al.27||EDI SD-OCT||Topcon, 3D-OCT 2000||50||•|
|Ikuno et al.28||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Li et al.29||EDI SD-OCT||Heidelberg, Spectralis||25||•|
|Li et al.30||EDI SD-OCT||Heidelberg, Spectralis||25||•|
|Li et al.9||EDI SD-OCT||Heidelberg, Spectralis||?||•|
|Lindner et al.31||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Margolis and Spaide32||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Matsuo et al.37||EDI SD-OCT||Heidelberg, Spectralis||100||•|
|Matsuo et al.37||EDI SD-OCT||Topcon, 3D-OCT 1000||50||•|
|Read et al.5||EDI SD-OCT||Optopol, Copernicus SOCT-HR||?||•|
|Shin et al.34||EDI SD-OCT||Topcon, 3D-OCT 2000||8||•|
|Sim et al.35||EDI SD-OCT||Heidelberg, Spectralis||?||•|
|Alonso-Caneiro et al.38||EDI SD-OCT||Heidelberg, Spectralis||30||•|
|Lee et al.36||EDI SD-OCT||Carl Zeiss Meditec, Cirrus HD-OCT||20||•|
|Tian et al.40||EDI SD-OCT||Heidelberg, Spectralis||?||•|
|Esmaeelpour et al.25||SS-OCT||3D 1060 nm-OCT||?||•|
|Esmaeelpour et al.26||SS-OCT||3D 1060 nm-OCT||?||•|
|Hirata et al.12||SS-OCT||Topcon, SS-OCT||?||•|
|Ikuno et al.28||SS-OCT||1060 nm-OCT||?||•|
|Matsuo et al.37||SS-OCT||Atlantis, DRI OCT-1||96||•|
|Ruiz-Moreno et al.33||SS-OCT||Topcon, SS-OCT||96||•|
|Usui et al.8||SS-OCT||Topcon, SS-OCT||32||•|
|Kajić et al.43||SS-OCT||3D 1060 nm-OCT||?||•|
|Zhang et al.42||SS-OCT||Topcon, SS-OCT||?||•|
|Torzicky et al.13||PS-OCT||1040 nm-OCT||?||•|
In general, manual delineation of the OCB, even with assistance of built-in software (calipers), must cope with serious difficulties as wide intra- and inter-observer variability and large time effort (cf. Ref. 28, p. 5539; Ref. 37, p. 7635). The crucial point, however, is the required quality of the underlying OCT data. Table 6 shows that, for manual OCB segmentation, highly averaged data are widely preferred. The results of the present study confirm these observations, showing a considerable variation in the five manual delineations. Thus, for moderately averaged SD-OCT data, a stable automated method as the proposed one seems to be clearly preferable to manual segmentation.
For an exemplary comparison with graph-theoretical or statistical automated segmentation methods, we choose the papers 38, 39, and 43 (the setting of Ref. 39 is closest to the present study). By overall inspection, our approach leads to results that are fairly well comparable with the aforementioned studies, even though Refs. 38 and 43 used data generated by special imaging modules, containing a considerably lower amount of noise. The errors in OCB segmentation and choroidal thickness measurement obtained in Refs. 38 and 39 fit into the same range as the ones obtained in the present study. A detailed quantitative comparison with the results from Ref. 43 is not possible due to the different error measure used there. Since their method involves model parameters generated from the examination of a training set, it is unclear whether the measurements are externally repeatable. Furthermore, for this and the method from Ref. 38, it may be questioned what constitutes the lower bound of image quality limiting their application. For graph-theoretical methods, a general challenge is to decide whether the possible occurrence of long-range correlations influences local statistics.
Comparison of manual and automated ILM/RPE boundary segmentation. In order to demonstrate its reliability, the quadratic measure filter method has been applied purposefully even to the “simple” detection of the ILM and RPE/choroid boundaries. In Table 7, the information from Table 2 is summarized, showing for every column the mean and the standard deviation of its 50 entries.
Summary: results of automated ILM and RPE detection. Quantities defined as in Table 2.
|ILM‾σ (μm)||ILM¯-|Δ|am (μm)||ILM¯-Δam (μm)||RPE‾σ (μm)||RPE¯-|Δ|am (μm)||RPE¯-Δam (μm)|
|Mean of the entries from Table 2||2.91||3.24||3.24||7.16|
|Std. dev. of the entries Table 2||(n/a)||(n/a)|
More precisely, we find that and are less than 1 px () for 45 of the 50 and 50 of the 50 investigated scans, respectively. Further, is less than 2 px for 34 of the 50 and less than 3 px for 45 of the 50 scans while is less than 2 px for 38 of the 50 and less than 3 px for 45 of the 50 scans. The accuracy of these results is comparable with the mean absolute errors mentioned in the literature, e.g., Ref. 18, p. 52, Fig. 22 C (ILM boundary: , RPE boundary: , automated detection versus mean of three manual delineations), Ref. 21, p. 1752, Table 1, third column (ILM boundary: , RPE boundary: , automated detection versus single manual delineation) and Ref. 38, p. 2809, Table 1 (RPE boundary: in two samples, automated detection versus single manual delineation). Our results are of the same magnitude for the ILM boundary and slightly larger for the RPE/choroid boundary. For the latter, the mean absolute error reduces to and the mean error to when the five outliers (Table 2, lines 9, 11, 15, 18, and 32) are excluded. We therefore summarize that the reliability of the quadratic measure filter method is well confirmed by our results.
Comparison of manual and automated OCB segmentation. As mentioned above, the manual delineation of OCB shows high variability between the five graders; see Table 3. The data sample contains a number of scans showing the OCB with a very weak contrast, leading to considerable disagreement between the manual segmentations. In the following Table 8, the information from Table 3 is summarized, showing the mean and the standard deviation of its 50 entries for columns 1 to 4.
Summary: results of automated OCB detection. Quantities defined as in Table 3.
|OCB‾σ (μm)||OCB¯-|Δ|am (μm)||OCB¯-Δam (μm)||Perc. (%)|
|Mean of the entries from Table 3||43.67||35.76||83.3|
|Std. dev. of the entries from Table 3||(n/a)|
Our automated OCB determination satisfies for 38 of the 50 investigated scans. The reliability of the method is shown as well by the fact that the inclusion is satisfied on more than 75% of the interval length for 41 of the 50 scans. In Ref. 39, analyzing data with a comparable amount of noise, a mean absolute error of and a mean error of were obtained (ibid., p. 1724, Table 2, automated detection versus consensus of two manual delineations). In Ref. 38, using EDI SD-OCT data with higher averaging, mean absolute errors of and and mean errors of and were documented (ibid., p. 2809, Table 1, two samples, automated detection versus single manual delineation). Our values, as summarized in Table 8, are slightly larger.
Comparison of manual and automated choroidal thickness measurements. The results from Tables 4 and 5 are summarized in the following Tables 9–11. More precisely, the first four columns of Table 9 contain the mean and the standard deviation of the 50 entries from Table 4, columns 6 to 9. In column 5, we appended mean and standard deviation of the absolute values of the entries from Table 4, column 9. Analogously, the first four columns of Tables 10 and 11 display the mean and the standard deviation of the 50 entries from Table 5, columns 1 to 4 and 5 to 8, respectively. In addition, we show the mean and standard deviation of the absolute values of the errors from Table 5, columns 4 and 8, respectively.
Summary: choroidal thickness, central foveal region. Quantities defined as in Table 4. Automated measurements in boldface.
|CHCm (μm)||CHCσ (μm)||CHCa (μm)||CHC-Δam (μm)|||CHC-Δam| (μm)|
|Mean of the entries from Table 4||283.56||55.93||279.29||25.64|
|Std. dev. of the entries from Table 4||(n/a)|
Summary: choroidal thickness, left outer 1 mm region. Quantities defined as in Table 5. Automated measurements in boldface.
|CHLm (μm)||CHLσ (μm)||CHLa (μm)||CHL-Δam (μm)|||CHL-Δam| (μm)|
|Mean of the entries from Table 5||272.83||50.71||275.65||26.76|
|Std. dev. of the entries from Table 5||(n/a)|
Summary: choroidal thickness, right outer 1 mm region. Quantities defined as in Table 5. Automated measurements in boldface.
|CHRm (μm)||CHRσ (μm)||CHRa (μm)||CHR-Δam (μm)|||CHR-Δam| (μm)|
|Mean of the entries from Table 5||248.58||40.34||254.28||19.88|
|Std. dev. of the entries from Table 5||(n/a)|
Again, the outcomes of the automated procedure fit excellently into the standard deviation intervals of the values generated from the manual delineations. For the central foveal region, we find for 46 of the 50 investigated scans; for the left 1-mm region, we have for 46 of the 50 scans, and for the right 1-mm region, we get for 44 of the 50 scans.
The errors observed in Tables 9–11 compete fairly well with Refs. 38 and 39, who averaged choroidal thickness measurements over the full length of the B-scans. Reference 38 obtained mean absolute errors of and and mean errors of and (ibid., p. 2810, Table 2, two samples, automated detection versus single manual delineation) while Ref. 39 found a mean error of (ibid., p. 1724, Table 2, automated detection versus consensus of two manual delineations). The agreement between manual and automated choroidal thickness measurements for the central foveal region is visualized by Bland–Altman plots;6566.–67 see Fig. 4. The shown examples are based on Table 4 (central foveal region). The agreement between the mean of the observers and the automated method is excellent. The five plots of automated versus single manual thickness measurements confirm considerable inter-observer variation when measuring the same subjects. Nevertheless, the automated method fits well into the 95% boundaries in all cases. Once more it becomes clear that, for data with the given image quality, the automated method is superior to a single manual determination.
Discussion of the obtained thickness values. Although not the focus of the present study, the obtained thickness measurements will be briefly set in position. The automatically derived choroidal thicknesses within our study are in line with previously obtained results from the literature; see Ref. 5, p. 3589, Table 1, column 1 (mean of eight foveal locations, , children of ages 4 to 12), Ref. 29 p. 552, Table 1 (central foveal region, , boys of ages 11 to 12 and , girls of ages 11 to 12), Ref. 30, p. 621, Table 3 (macular region, , , and , three groups of children of ages 11 to 12) and, Ref. 33, p. 353 (macular region, , children of ages 3 to 18). Exemplarily, we add a plot displaying the relation between the choroidal thicknesses for the central foveal region and the axial length (Fig. 5), which shows a clustering of boys and girls analogously to Ref. 29, p. 553, Fig. 1. Due to the small number of probands, we refrain from claiming other trends from this illustration. In a forthcoming study, we plan to analyze a larger sample from the LIFE cohort study.
Limitations of the Proposed Method; Possible Errors
Methodological limitations and possible errors. Within the proposed algorithm, there are at least two possible sources of systematic errors. First, choroidal thickness is measured in the vertical direction within the scan instead of the (geometrically) normal direction to the obtained RPE boundary, thus neglecting the curvature of the layer boundary. Second, Eq. (7) involved in the determination of the left and right 1-mm measurement regions is based on a model eye of an adult instead of a child.
Automated layer detection may produce errors, which a human observer would avoid. In rare cases, the feature recognition routine involved in the detection of the RPE boundary misidentifies parts of the next hyperreflective layer above the RPE for the RPE boundary itself, thus overestimating the resulting choroidal thickness (this error accounts for the outliers of and in Table 2, lines 9, 11, 15, 18, and 32). If druses were present (not occurring in the studied sample), similar errors would be possible.
The line averaging of 10 used in the study constitutes a lower bound for the applicability of the method, which becomes infeasible for data corrupted with stronger noise. It is to be expected that for older adult subjects, more averaging is required due to naturally occurring opacifications in the visual pathway.
Limitations due to choroid physiology. As mentioned in the introduction, choroidal thickness can be short-term modulated, even following a diurnal rhythm. As the prevailing part of the literature, the present study is unable to discern such effects when performing single OCT measurements.
Observed structure of the OCB compared with histology. In difference to Refs. 36 and 38 as well as to most studies performing manual OCB delineation, we observed a relatively bumpy structure of the OCB in most of our automated segmentations [the result shown in Fig. 1(b) is typical]. In the following, we will shortly discuss whether our observations correspond to the histological composition of the suprachoroid and choroid/sclera border.
Note first that, in the case of manual segmentation, the OCB is usually fixed by marking a comparatively small number of measurement points [see, e.g., Ref. 4, p. 326, Fig. 1 (11 points), Ref. 24, p. 9556, Fig. 1 (9 points), Ref. 28, p. 5537, Fig. 1 (6 points), and Ref. 31, p. 876 (25 points per volume scan)] and possible subsequent spline interpolation. Necessarily, boundary lines resp. thickness maps resulting from such procedures appear very smooth. In contrast, there are several studies with automated segmentation reporting a bumpy structure of the OCB [we mention Ref. 13, p. 7570, Fig. 2, Ref. 42, p. 3204, Fig. 4(d), and Ref. 43, p. 99, Fig. 5(b)].
Histological investigation of the human suprachoroid shows an irregular distribution of features with low and high reflectance in visible light (melanocytes resp. translucent parts and collagen fibers) throughout this transition layer [Fig. 6(a)]. Cross-section reveals a sharply defined histological border between the suprachoroid and the collagen bundles of the sclera which, however, does not necessarily correspond to an optically distinguishable brightness jump [absent in Fig. 6(b), present in Fig. 6(c)]. Therefore, it is likely that the histological border may not be represented as a separate reflectivity feature within the OCT data. Consequently, our method provides a highly reliable inner approximation of the histologically determined choroid boundary but possibly underestimates the choroidal thickness. If so, the same would be true for the vast majority of thickness measurements presented in the literature. In view of Fig. 6(a), the bumpy appearance of the detected OCB may well correspond to the scattered distribution of suprachoroidal tissue regions with sufficiently high concentration of melanocytes.
A novel method for the automated detection of the OCB within SD-OCT image data based on an image model within the space of BV functions and the application of quadratic measure filters was established and proved to be reliable. Additionally, the method was successfully applied to the automated segmentation of the ILM and RPE/choroid boundaries. A particular advantage of the presented approach is its capability of dealing with data generated without special imaging modules and moderate averaging only. For such data, as available in the LIFE Child and LIFE Adult cohort studies, the method is clearly preferable to manual OCB segmentations and subsequent thickness derivations. The presented algorithms can be externally repeated independent of the sample used in the study.
Edge Detection within Noisy Grayscale Images
Our approach fits into a mathematical framework where a grayscale image containing the reflectivity information is understood as a function rather than as a “discrete” entity. Then the acquisition of OCT data can be modeled as68, p. 1208. Of course, the edge structure in is distorted by the noise term as well. Consequently, the recognition of edges requires a partial removal of the noise at least.
In large parts of the literature, denoising is performed first by subjecting the data to a suitable diffusion process (see Ref. 69, pp. 94 ff., for a closer mathematical analysis of diffusion-based image restoration), which is followed by separate, more or less sophisticated edge detection steps. In another largely pursued approach, the formal reconstruction of , which is unstable due to the presence of the noise term, is replaced by a regularized variational or optimal control problem in suitable function spaces. Often, the cost functionals are designed particularly for the suppression of speckle noise.7071.–72 While solving such problems, a simultaneous edge detection is often possible, e.g., by use of Ambrosio–Tortorelli type functionals (see Ref. 69, p. 169 f.) or by incorporation of gradient constraints.73
Functions of Bounded Variation
In the present study, we model an OCT image as an integrable function of bounded variation (BV) defined on a rectangle . A BV function must satisfy the condition74, pp. 166 ff., for a rigorous mathematical treatment). The use of this function space in image processing has a long and successful tradition (see Ref. 75, pp. 174 ff.) since, on one hand, a notion of weak derivatives is available, and on the other hand, functions of BVs are allowed to have jumps. More precisely, the edge structure within the image is defined as the union of finitely many piecewise smooth curves, along which the function admits jumps of the height , . In the following, we will outline how this edge structure can be identified by the approach of quadratic measure filters.45,4748.–49 From the methods described above, this approach differs in the fact that the smoothing kernel applied to the noisy data forms an integral part of the edge detector itself.
Edge Detection by Quadratic Measure Filters
In mathematical terms, the weak partial derivative of a BV function can be identified with the sum of three signed measures , , and : the first one (which is absolutely continuous with respect to the 2-D Lebesgue measure ) representing the usual partial derivative where it exists but measuring the jump set with zero, the second one encoding the position and height of possible jumps while the third part (the so-called Cantor part) can be neglected in medical image processing. Consequently, in order to identify the edge structure within , one has to access the measure . Indeed, this is possible by the following theorem yielding an (implicit) description of in terms of a limit relation for measures.
(Ref. 48, p. 2, Theorem 1; cf. also Ref. 45, p. 701 f., Theorem 3.1.) We choose a nonnegative, symmetric kernel , which is concentrated on the unit ball and normed by . With the aid of , we build a family of smoothing kernels , . Then for every integrable function of bounded variation, the following limit relation holds true:
Since the limit in Eq. (33) is formed with respect to (weak*-) convergence in the space of signed measures, the assertion takes the form of a variational equation, which holds true for a whole set of test functions . The integral on the left-hand side of Eq. (33) is termed a “quadratic measure filter” and can be understood as the approximation of a quantity, which is concentrated on the jump set of and measures the square of the jump height. Note that the theorem remains true for the Gauss kernel and noncompactly supported test functions as well, cf. Ref. 47.
Assuming that an admissible test function is supported on a subset , Eq. (33) implies the alternative2.3 are derived by insertion of a fixed test function supported on and subsequent discretization of the integral (here the function is replaced by the pixeled data ). Then, instead of passing to the limit, one may inspect the local maxima of the resulting detectors for fixed while moving the test function and its support through appropriate subsets of . Namely, the maxima search will be constrained to the intersection of the sets or with narrow vertical stripes, thus accounting for the mainly horizontal orientation of the layer structure within OCT data.
This publication is supported by LIFE Leipzig Research Center for Civilization Diseases, Leipzig University. LIFE is funded by means of the European Union, the European Social Fund (ESF), the European Regional Development Fund (ERDF) and the Free State of Saxony within the framework of the excellence initiative. The Integrated Research and Treatment Center Adiposity Diseases is funded by the German Federal Ministry of Education and Research (Grant No. 01EO1501). Marcus Wagner gratefully acknowledges support by the Max-Planck-Gesellschaft for granting a stay at the MPI for Mathematics in the Sciences, Leipzig, in 2015 and 2016. The authors would like to express their sincere thanks to Heike Lange (Leipzig University Hospital, Department of Ophthalmology/LIFE Leipzig Research Center for Civilization Diseases, Leipzig University) for carrying out the presented optical coherence tomography measurements as well as obtaining the axial length and corneal radii values employed and to Claudia Holzhey (Heidelberg Engineering GmbH, Heidelberg) for her support with the implementation of the measurement technique as part of the LIFE Child study. Furthermore, the authors would like to thank the following researchers for carrying out manual segmentations for this paper: Anna Xenia Bestehorn (Münster), Sophia Scheibe (Leipzig), and Anja Schirmer (Berlin).
Patrick Scheibe received his MSc degree in computer science in 2007 from Leipzig University. For the following 9 years, he was head of the image processing core unit at the Translational Centre for Regenerative Medicine Leipzig. Since 2016, he is leading data analyst at the Saxonian Incubator for Clinical Translation. His research interests focus on the quantification and analysis of biomedical experiments. Currently, he concentrates on modeling of foveal structures based on OCT measurements.