CuBe : parametric modeling of 3 D foveal shape using cubic Bézier

Optical coherence tomography (OCT) allows three-dimensional (3D) imaging of the retina, and is commonly used for assessing pathological changes of fovea and macula in many diseases. Many neuroinflammatory conditions are known to cause modifications to the fovea shape. In this paper, we propose a method for parametric modeling of the foveal shape. Our method exploits invariant features of the macula from OCT data and applies a cubic Bézier polynomial along with a least square optimization to produce a best fit parametric model of the fovea. Additionally, we provide several parameters of the foveal shape based on the proposed 3D parametric modeling. Our quantitative and visual results show that the proposed model is not only able to reconstruct important features from the foveal shape, but also produces less error compared to the state-of-the-art methods. Finally, we apply the model in a comparison of healthy control eyes and eyes from patients with neuroinflammatory central nervous system disorders and optic neuritis, and show that several derived model parameters show significant differences between the two groups. © 2017 Optical Society of America OCIS codes: (100.0100) Image processing; (100.2960) Image analysis; (110.4500) Optical coherence tomography; (000.4430) Numerical approximation and analysis. References and links 1. D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, C. Puliafito, and a. et, “Optical coherence tomography,” Science 254, 1178–1181 (1991). 2. N. M. Roth, S. Saidha, H. Zimmermann, A. U. Brandt, J. Isensee, A. Benkhellouf-Rutkowska, M. Dornauer, A. A. Kühn, T. Müller, P. A. Calabresi, and F. Paul, “Photoreceptor layer thinning in idiopathic parkinson’s disease,” Movement Disorders 29, 1163–1170 (2014). 3. S. Stricker, T. Oberwahrenbrock, H. Zimmermann, J. Schroeter, M. Endres, A. U. Brandt, and F. Paul, “Temporal retinal nerve fiber loss in patients with spinocerebellar ataxia type 1,” PLoS ONE 6, e23024 (2011). 4. T. Oberwahrenbrock, M. Ringelstein, S. Jentschke, K. Deuschle, K. Klumbies, J. Bellmann-Strobl, J. Harmel, K. Ruprecht, S. Schippling, H.-P. Hartung, O. Aktas, A. U. Brandt, and F. Paul, “Retinal ganglion cell and inner plexiform layer thinning in clinically isolated syndrome,” Multiple Sclerosis Journal 19, 1887–1895 (2013). 5. L. J. Balk, P. Tewarie, J. Killestein, C. H. Polman, B. Uitdehaag, and A. Petzold, “Disease course heterogeneity and oct in multiple sclerosis.” Multiple sclerosis (Houndmills, Basingstoke, England) 20, 1198–1206 (2014). 6. F. Pache, H. Zimmermann, J. Mikolajczak, S. Schumacher, A. Lacheta, F. C. Oertel, J. Bellmann-Strobl, S. Jarius, B. Wildemann, M. Reindl, A. Waldman, K. Soelberg, N. Asgari, M. Ringelstein, O. Aktas, N. Gross, M. Buttmann, T. Ach, K. Ruprecht, F. Paul, A. U. Brandt, and in cooperation with the Neuromyelitis Optica Study Group (NEMOS), “Mog-igg in nmo and related disorders: a multicenter study of 50 patients. part 4: Afferent visual system damage after optic neuritis in mog-igg-seropositive versus aqp4-igg-seropositive patients.” J. Neuroinflammation 13, 282 (2016). 7. E. Schneider, H. Zimmermann, T. Oberwahrenbrock, F. Kaufhold, E. M. Kadas, A. Petzold, F. Bilger, N. Borisow, S. Jarius, B. Wildemann, K. Ruprecht, A. U. Brandt, and F. Paul, “Optical coherence tomography reveals distinct patterns of retinal damage in neuromyelitis optica and multiple sclerosis,” PLOS ONE 8, 1–10 (2013). 8. F. C. Oertel, J. Kuchling, H. Zimmermann, C. Chien, F. Schmidt, B. Knier, J. Bellmann-Strobl, T. Korn, M. Scheel, A. Klistorner, K. Ruprecht, F. Paul, and A. U. Brandt, “Microstructural visual system changes in Vol. 8, No. 9 | 1 Sep 2017 | BIOMEDICAL OPTICS EXPRESS 4181 #295770 Journal © 2017 https://doi.org/10.1364/BOE.8.004181 Received 15 May 2017; revised 26 Jul 2017; accepted 27 Jul 2017; published 22 Aug 2017 AQP4-antibody–seropositive NMOSD,” Neurology Neuroimmunology Neuroinflammation 4, e334 (2017). 9. T. Oberwahrenbrock, M.Weinhold, J. Mikolajczak, H. Zimmermann, F. Paul, I. Beckers, and A. U. Brandt, “Reliability of intra-retinal layer thickness estimates,” PloS one 10, e0137316 (2015). 10. R. Kafieh, H. Rabbani, and S. Kermani, “A review of algorithms for segmentation of optical coherence tomography from retina.” J. Medical Signals and Sensors 3, 45–60 (2013). 11. S. Niu, Q. Chen, L. de Sisternes, D. L. Rubin, W. Zhang, and Q. Liu, “Automated retinal layers segmentation in sd-oct images using dual-gradient and spatial correlation smoothness constraint,” Computers in Biology and Medicine 54, 116–128 (2014). 12. Y. M. Cha and J. H. Han, “High-accuracy retinal layer segmentation for optical coherence tomography using tracking kernels based on gaussian mixture model,” IEEE J. Sel. Top. Quantum Electron. 20, 32–41 (2014). 13. L. Fang, D. Cunefare, C. Wang, R. H. Guymer, S. Li, and S. Farsiu, “Automatic segmentation of nine retinal layer boundaries in oct images of non-exudative amd patients using deep learning and graph search,” Biomed. Opt. Express 8, 2732–2744 (2017). 14. S. Tick, F. Rossant, I. Ghorbel, A. Gaudric, J.-A. Sahel, P. Chaumet-Riffaud, and M. Paques, “Foveal shape and structure in a normal population,” Investigative Ophthalmology and Visual Science 52, 5105 (2011). 15. T. Y. P. Chui, Z. Zhong, H. Song, and S. A. Burns, “Foveal avascular zone and its relationship to foveal pit shape,” Optometry and Vision Science 89, 602–610 (2012). 16. D. X. Hammer, N. V. Iftimia, R. D. Ferguson, C. E. Bigelow, T. E. Ustun, A. M. Barnaby, and A. B. Fulton, “Foveal fine structure in retinopathy of prematurity: An adaptive optics fourier domain optical coherence tomography study,” Investigative Ophthalmology and Visual Science 49, 2061 (2008). 17. B. K. McCafferty, M. A. Wilk, J. T. McAllister, K. E. Stepien, A. M. Dubis, M. H. Brilliant, J. L. Anderson, J. Carroll, and C. G. Summers, “Clinical insights into foveal morphology in albinism,” J. Pediatric Ophthalmology Strabismus 52, 167–172 (2015). 18. Y. Barak, M. P. Sherman, and S. Schaal, “Mathematical Analysis of Specific Anatomic Foveal Configurations Predisposing to the Formation of Macular Holes,” Investigative Ophthalmology Visual Science 52, 8266–8270 (2011). 19. B. Nesmith, A. Gupta, T. Strange, Y. Schaal, and S. Schaal, “Mathematical Analysis of the Normal Anatomy of the Aging FoveaAnatomy of the Aging Fovea,” Investigative Ophthalmology & Visual Science 55, 5962–5966 (2014). 20. A. M. Dubis, J. T. McAllister, and J. Carroll, “Reconstructing foveal pit morphology from optical coherence tomography imaging,” Br. J. Ophthalmol. 93, 1223–1227 (2009). 21. P. Scheibe, M. T. Zocher, M. Francke, and F. G. Rauscher, “Analysis of foveal characteristics and their asymmetries in the normal population,” Experimental Eye Research 148, 1–11 (2016). 22. P. Scheibe, A. Lazareva, U.-D. Braumann, A. Reichenbach, P. Wiedemann, M. Francke, and F. G. Rauscher, “Parametric model for the 3d reconstruction of individual fovea shape from OCT data,” Experimental Eye Research 119, 19–26 (2014). 23. M. A. Wilk, J. T. McAllister, R. F. Cooper, A. M. Dubis, T. N. Patitucci, P. Summerfelt, J. L. Anderson, K. E. Stepien, D. M. Costakos, T. B. Connor, Jr, W. J. Wirostko, P.-W. Chiang, A. Dubra, C. A. Curcio, M. H. Brilliant, C. G. Summers, and J. Carroll, “Relationship between foveal cone specialization and pit morphology in albinismfoveal cone specialization and pit morphology,” Investigative Ophthalmology and Visual Science 55, 4186 (2014). 24. L. Liu, W. Marsh-Tootle, E. N. Harb, W. Hou, Q. Zhang, H. A. Anderson, T. T. Norton, K. K. Weise, J. E. Gwiazda, L. Hyman, and the COMET Group, “A sloped piecemeal Gaussian model for characterising foveal pit shape,” Ophthalmic and Physiological Optics 36, 615–631 (2016). 25. Y. Ding, B. Spund, S. Glazman, E. M. Shrier, S. Miri, I. Selesnick, and I. Bodis-Wollner, “Application of an oct data-based mathematical model of the foveal pit in parkinson disease,” Journal of Neural Transmission 121, 1367–1376 (2014). 26. G. E. Farin, Curves and Surfaces for Computer-Aided Geometric Design: A Practical Code (Academic Press, Inc., Orlando, FL, USA, 1996), 4th ed. 27. R. Ihaka and R. Gentleman, “R: A language for data analysis and graphics,” Journal of Computational and Graphical Statistics 5, 299–314 (1996). 28. F. C. Oertel, H. Zimmermann, J. Mikolajczak, M. Weinhold, E. M. Kadas, T. Oberwahrenbrock, F. Pache, J. BellmannStrobl, K. Ruprecht, F. Paul, and A. U. Brandt, “Contribution of blood vessels to retinal nerve fiber layer thickness in nmosd,” Neurology Neuroimmunology Neuroinflammation 4 (2017).


Introduction
Optical coherence tomography (OCT) is a non-invasive imaging modality based on low-coherence interferometry [1].OCT has become an essential diagnostic tool in ophthalmology for its capability to image the retina at micrometer resolution.
Recently, the scope of retinal OCT imaging has extended towards neurologic diseases.As part of the central nervous system, the retina comprises a similar cellular composition as the brain, and many neurologic disorders thus affect the retina.OCT is able to detect retinal changes in neurodegenerative disorders like Parkinson's disease [2] and cerebellar ataxias [3] as well as autoimmune neuroinflammatory disorders like clinically isolated syndrome [4], multiple sclerosis [5], MOG-ab-positive encephalomyelitis [6], and neuromyelitis optica spectrum disorders [7,8].
Traditionally, OCT has been used to detect and describe macroscopic changes.However, neurologic diseases often lead to only small retinal changes, which rather need to be measured than can be seen directly.In the following, we will use the term quantitative OCT for this approach.Commonly, quantitative OCT determines the size of retinal structures as thickness or volume in a defined area of interest [9].Structures of interest for measurements are the optic nerve head and the macula (Fig. 1).For example, one of the quantitative OCT parameters is the peripapillary retinal nerve fiber layer thickness (pRNFL), which is usually measured in a 12°ring scan around the optic nerve head.Quantitative measurements of macular changes were initially defined by, for example, total macular thickness (TMV) or total retinal thickness (TRT), which measure the full retinal thickness across the macula [9].Recently, intra-retinal layer segmentation at the macula has been developed to allow quantification of intraretinal thickness or volume changes, for example, of the ganglion cell layer (GCL) [10][11][12][13].
Foveal shape analysis is a promising approach for quantitative OCT next to thickness or volume measurements.The few studies in this regard can be divided into two categories: data-driven, which directly employ segmentation lines from the OCT scan itself to compute various metrics, and model-driven, which use mathematical constructs to create a representation of the macula and fovea in order to calculate different features directly from the model.

Data-driven approaches
All data-driven approaches known to the authors are based on the analysis of 2D images (B-scans) with results computed from one B-scan or averaged over two B-scans.
To analyze the variability of the healthy foveal shape and to investigate the relationship between this structure and the foveal avascular zone (FAZ), Tick et al. computed several parameters: pit depth, central foveal thickness, maximal retinal thickness, pit diameter, pit cross-sectional area, and the foveal inner retinal area.Using these parameters, the authors could show that the healthy foveal structure strongly correlates with its neurovascular structure [14].
Similarly, Chiu et al. introduced the parameters foveal photoreceptor thickness and foveal width, and found a large individual variation of FAZ size and shape in healthy retinae and a negative correlation between FAZ diameter and foveal thickness [15].
Other studies investigated the way retinal structures change in the context of prematurity [16] or albinism [17] using fovea shape parameters like pit depth.

Model-driven approaches
One of the first attempts to mathematically describe the foveal shape was published by Barak et al. [18].Their approach employs an automated symbolic regression software that fits a section of the foveal profile around the center which was able to detect different patterns in the premacular hole foveal configurations and normal foveal configurations.A similar approach was used by Nesmith et al. to characterize changes that occur in the foveal anatomy with aging in a large number of OCT scans [19].
The first more general mathematical model was created by Dubis et al., who used Difference of Gaussians (DoG) function for a 2D fit of several radial scans passing through the lowest point of a macular scan [20].The method derived three pit metrics (diameter, depth and slope) from a fitted model, which was symmetrical with respect to the fovea center.The model's symmetry was problematic, because of the fovea's asymmetry, which stems mostly from differences between the nasal and temporal retinal nerve fibers [14,21].The model was thus not accurate enough for modeling the asymmetric nature of the fovea and had problems covering the wide range of fovea shapes in the healthy population, due to other model fitting constraints, which were e.g.not able to model an extended flat part at the foveal pit in many healthy foveae [14,21].
Scheibe et al. applied a different approach to overcome these drawbacks by creating a more flexible 2D model, which is fitted to each supporting direction in a circular region around the foveal center [22].This method used a Gaussian like basis function (exp(−x γ )), when γ = 2, it would be a Gaussian basis function otherwise it would be an exponential function.The parameter γ defines the shape of the exponential basis function.The authors derived five parameters analytically: mean retinal thickness inside a radius of 1mm, foveal bowl area, retinal radius, and maximum height of the foveal rim.By applying this model to a large cohort of healthy subjects, the authors were able to confirm and further investigate the asymmetric nature of the fovea.Although, capable of modeling a large variety of data, and having a smaller root mean square error (RMSE) than previous approaches, the method has difficulties deriving 3D parameters such as volume due to the complexity of the computation because of the parameter γ.The RMSE is a standard metric to measure the model error.The low RMSE value indicates the high fidelity of a model.
Wilk et al. [23] used a similar algorithm described in [20] and applied it on a volume scan by extracting 180 radially oriented slices through the foveal center.In addition to the parameters defined in [20], the foveal pit volume was computed by calculating the space between the internal limiting membrane surface and the top of the foveal pit.
Liu et al. created a different 2D model, which uses a sloped piecemeal Gaussian function (SPG) to model the asymmetry and the foveal flatness [24].The authors tested their method on a large number of macular scans and showed good RMSE compared to other 2D methods.The main drawback consists in the use of only two scans for each subject, a vertical and a horizontal one, which limits information on only these sections and prohibits generation of 3D metrics like volume.
Ding et al. developed an approach, which intends to model the surface of inner retinal layers using a 3D model-fit based DoG and a second order polynomial [25].The authors used the model coefficients to discriminate between Parkinson patients and a control group.As a drawback, the model coefficients' relation to foveal shape itself is difficult to interpret.Additionally, although the model takes into account the difference in the slope for horizontal and vertical directions, it does not model the asymmetry in the anatomical ones (temporal-nasal), nor does it capture fovea shapes with very flat pits accurately.
Objective Against this background, the goal of this study was to develop a robust, modeldriven 3D macula shape analysis method, which can be computed from standard macular volume OCT scans and from which foveal and macular shape metrics can be derived.In addition to available macular thickness measurements, our approach allows a detailed analysis of foveal shape, including depth, diameter, slope, area and volume of different regions, as well as pit shape analysis.

Method theory
Most of the model driven state of art methods assume that the foveal pit shape is quite similar to the Gaussian function.Therefore, these methods use the Gaussian basis function for the mathematical modeling of the foveal pit [20,22,24,25].However, this assumption is not entirely true as we can see in Fig. 2 that the pit shape varies a lot from the Gaussian shape.Modeling the flatness of foveal pit is difficult with Gaussian basis function.The two regions of the Gaussian function (around and away from mean value) are dependent on each other and on the standard deviation of the given Gaussian function.However, the shape of the interior (pit area) and exterior (rim area) regions of fovea are independent of each other.For example, similar kinds of foveal pits can have significantly different rim heights.In summary, we can say that the Gaussian function is not an optimal basis function to describe macular and foveal shape.
To overcome these problems, we introduce a cubic Bézier based robust and flexible basis function, which is able to encounter all possible variations of the fovea shape.We fit the cubic Bézier polynomial in interior and exterior regions independently using a least square optimization with invariant features of fovea which will be discussed in the coming sections.

Invariant features of the fovea
Foveal shape in healthy population varies [21].Figure 2, shows a selection of different foveal pit shapes from our cohort to exemplify this.In order to reconstruct all these variations, some invariant features of the shape are needed.Tangents at critical points of macula represent reliable and stable features.The critical points refer to the lowest point of the pit and the highest point around this region as shown in Fig. 3.The slope at these critical points is always zero and macula will have horizontal tangents at these points.

Cubic Bézier
Bézier curves were introduced by Dr. Pierre Bézier in early 1960s.A parametric Bézier curve of polynomial degree n is defined as: where the P i,n are the control points and B i,n are the Bernstein polynomials.
For a Bézier curve of degree n, there are n + 1 control points.From Eq. ( 1), one can see that the Bézier curve is a weighted average of control points where weights are defined using Bernstein polynomials.For n = 3, Eq. ( 1) becomes the cubic Bézier equation: (3) Cubic Bézier curve will have four control points and is tangent to the first and last control points P 0 and P 3 respectively [26].The relationship between control points, in terms of distance α, β and unit tangent directions T 01 and T 23 can be written as: Distance parameters α and β represent the distances between end control points (P 0 , P 3 ) and inner control points (P 1 , P 2 ) respectively.In our curve fitting algorithm, parameters α and β of each Bézier segment are used as shape parameters for an optimal curve fitting.The direction of a Bézier curve at its endpoints is uniquely determined by the tangent vector.Thus, by choosing the same tangent vector for two adjacent Bezier segments tangent continuity (geometrical continuity G 1 ) at each junction is assured, and thus throughout the whole composite spline curve.

Least squares optimization
A least squares fitting approach is applied for each B-scan/radial scan using a cubic Bézier curve and tangents at critical points, which are derived from the invariant features of the macular shape.Let us consider c(x) as a central B-scan, which contains the minimal retinal thickness point as shown in Fig. 3(a).The central B-scan can then be decomposed into two interior (c 1 (x), c 2 (x)) and two exterior (c 0 (x), c 3 (x)) segments, where the right half of the central B-scan is defined as the segment between the points (x m , c(x m )) and (x e , c(x e )) as shown in Fig. 3(a).The maximum point of the right half of the central B-scan can be computed, as represented by (x 2 , c(x 2 )).Two segments can be made using the critical point (x 2 , c(x 2 )) and each segment was denoted as c I : [x m , x 2 ] → R for the interior, and c E : [x 2 , x e ] → R for the exterior segment as shown in Fig. 3(b).
For interior segment (c I ), end points, P 0 and P 3 , have horizontal tangent lines T = [1, 0] and inner control points, P 1 and P 2 , which are lying opposite to each other w.r.t.their tangent direction, as shown in Fig. 3(b).Thus, Eq. ( 4) can be modified to give: For the exterior segment (c E ), only one of the end points has a horizontal tangent and one of the inner control points (P 1 ) lies in the same horizontal tangent direction.Equation ( 4) for exterior segments then becomes: The other inner control point, P 2 , for the exterior segment does not have a horizontal tangent direction (T = [1, 0]); therefore the value of P 2 is optimized without using any invariant feature of the macula.In this proposed method, control points P i are defined in R 2 space.Next, the optimized values of α and β are computed to achieve the best-fit cubic Bézier for each segment of the scan.The optimization process begins from the interior segment.By using Eq. ( 3) and 5, the cubic Bézier equation for interior segment (c I ) can be modified as follows: If the corresponding segment has m data points and is represented by D I then the least squares energy function can be expressed as: where t i represents discrete values of t corresponding to the data points D I i and it is computed using the uniform parametrization (t i = i/m ∀i ∈ [0, 1, ..., m]).Equation ( 8) shows a quadratic energy function in α and β for a given parametric value of t and higher order polynomial in t for a given α and β.The minimization of the above energy function then becomes a multispace optimization, since there are three different variables to be minimized.The multi-space optimization procedure for t, α and β occurs using the following steps: 1. First, Eq. ( 8) is a quadratic energy function in α and β for a given t, so it can be solved by using the linear system of equations.An initial value is given to t ∈ [0, 1/m, 2/m, • • • , 1] using the uniform parametrization and corresponding α and β values are calculated by computing the first derivative of the given energy function w.r.t.α and β.
To compute the gradient of the given energy in terms of α: By using Eq. ( 9), we can conclude by stating the following relation: where the left hand side of Eq. (11) shows only α as an unknown variable, since the terms of β are omitted due to differential operations.However, β can also be computed using Eq.( 9) in the same manner.
2. Second, the optimized value of t is calculated using the gradient descent method with the values for α and β that were computed in the previous step.
Similar to Eq. ( 10), the derivative of the given energy function w.r.t t is computed as follows: where Q I (t i , α, β) represents the first derivative of cubic Bézier in terms of t.The optimized value of t is computed using the following equation: where s represents the step size of the gradient descent method and ti is the optimized value of t.In the next iteration, ti will be used as the initial value of t to compute the new α and β.Throughout the whole experimentation the step size was fixed to s = 0.1.
These two steps are iterated until the stable minimum of the energy function in Eq. ( 8) is obtained.Equation ( 13) depends not only on t but also on α, β, thus allowing stable minimum to be reached.The minimum of the given energy can be achieved in 400 iterations.After 400 iterations, there is no significant change in the optimization result.
For the exterior segment (c E ), our optimization procedure is different compared to the interior segments as only one end point has a horizontal tangent.Using Eq. ( 6) and Eq. ( 3), we can have a modified equation for optimization: Similar to the energy function in Eq. ( 8), here we have an energy function, which depends on t, P 2 and α and which is defined as: where D E represents data points corresponding to the external segment.The minimization of this given function is performed in a similar fashion as the interior segment parameter optimization.First, the derivative of the energy function w.r.t.P 2 and α, for a given initial uniform parameterized t, is computed.Then, the optimized t with newly obtained values for P 2 and α is determined.

Materials and methods
To evaluate the proposed method macular volume scan (25 All scans underwent quality control by an experienced rater.Automatic layer segmentation was performed with the device's software (Eye Explorer 1.9.10.0 with viewing module 6.0.9.0).The institute's imaging database only contains images derived from local studies that were approved by the local ethics committee at the Charité -Universitätsmedizin Berlin and were conducted following the Declaration of Helsinki in its currently applicable version.All computations were carried out using MATLAB 2016b (MathWorks, Inc., Natick, MA, USA).The statistical analysis was conducted using R version 3.3.2[27].ICC (inter class correlation) and GEE (Generalized estimating equation) were used for statistical measurements.The ICC measures the repeatability of the proposed parameters and GEE (p-value) shows the significance between the two groups of data for each of the proposed parameters.

Method pipeline
Figure 4 shows the pipeline of our algorithm.In order to import data into MATLAB, Heidelberg Spectralis OCT raw data format was exported from the device.This data contains additional to the image information, the coordinates of the inner limiting membrane (ILM) and the lower boundary of the retinal pigment epithelium layer denoted throughout our paper for simplicity reasons as RPE.The whole algorithm is implemented in the following steps: 1. ILM-RPE Computation and Minimal Foveal Point Detection: In the first step, the height difference between the ILM and the RPE of each volume scan is extracted.This represents the macular thickness surface.Using this difference has the advantage of removing the slant of the scan created at the measurement and/or by the curved shape of the eye.Let us consider a volume scan and the corresponding thickness profile represented as the graph function M : (x, y) → R, where (x, y) ∈ Ω and Ω represents our region of interest.A volume scan is the combination of A-scans and B-scans obtained from the OCT scanner.We assume that x and y represent A-scans and B-scans directions respectively.
To determine the fovea's center, a region Ω of 1 mm radius is taken from the surface centered at the fovea automatically detected by the OCT device.The information about this center point is included in the raw data export.In order to detect the lowest point of foveal surface, we look at the minima of this region.
where x m , y m are the coordinate of the minimum value M m of the volume scan.If several minimas are detected, then the median point of them is taken as the center of foveal pit.
2. Volume to Radial Sampling: This is the second step of the proposed method as shown in Fig. 4. For 3D shape analysis of fovea, information from the whole volume is needed, and therefore the scan is re-sampled into a radial one.The radial scans capture the foveal pit shape accurately as they have more samples near the center compared to the outer region.
Sampling from volume M to radial M p can be done in the following steps: (a) Create a polar grid M p (r, θ), centered at (x m , y m ) with zero height value.Radius and angle between the radial lines will be defined by the user.
(b) Compute the height value of M p (r, θ) using bilinear interpolation between the nearest four points of the M(x, y) which are closest to the corresponding (r, θ).
(c) Now, there are n = (2π/θ) radial scans represented as: These radial scans approximate the original volume scan as shown in Fig. 5(a).
In our experimentation, we choose r = 2mm and θ = 15 • for radial sampling.4. Cubic Bézier Fitting using Least Square Optimization: For each of the segments of the radial scan, a cubic Bézier with least square optimization as explained in section 2.3 is fitted.For the interior segment, Eq. ( 8), ( 9) and ( 12) is used to compute optimized α, β and t.Then, the interior segment using the optimized parameters is reconstructed by assigning these to Eq. ( 7).Similarly, for exterior segment, Eq. ( 16) is used and optimized α, P 2 are computed.Then the exterior segment is reconstructed by using Eq. ( 15).To get a complete 3D parameterized modeling of the fovea, a cubic Bézier has to be fitted to each of the re-sampled radial scans.

Parameters
In this section, we present several shape parameters for a volume scan using a cubic Bézier parameterization.Rim point of a radial scan is defined as the maximum height point in the corresponding radial scan.Let us consider that there are n number of radial scans re-sampled from a volume scan and (p 1 , p 2 , • • • , p n ) ∈ R 3 are the corresponding rim points.Figure 6 shows re-sampled radial scans from a volume scan divided into two parts: c(r, θ i ) and c(r, θ i + π).The notation used to denote the corresponding least square optimization fitted radial scan parts is A visual representation of some of the 3D parameters is shown in Fig. 8 and 9. Figure 7 shows two parameterized radial scans (Q(t, θ i ) and Q(t, θ i + π)) of a volume scan which differ by the angle π.Several basic parameters are defined based on these two radial scans.Q I and Q E represent parameterized interior and exterior segments respectively.Most of the 3D parameters are computed by first computing the values for each radial scans and then taking the average over these scans.Central Foveal Thickness (h c f t ): refers to the central foveal thickness which is defined as the minimum height of fovea at the center of the pit.From Eq. ( 17), h c f t can be written as: The h c f t of each radial scan is the same because (x m , y m ) is the center of radial and represents the beginning as well as the lowest point for each radial scan.Average Rim Height (h r ): is defined as average of maximum height in each radial scan of a volume (as shown in Fig. 3).Average rim height is written as:

Maximum
Rim Disk Area (A r ): To compute the rim disk area, the normal to the disk plane is calculated as shown in Fig. 8.The covariance analysis of all the rim points provides not only the disk normal but also shape information as described below.Rim points covariance matrix can be computed as: Rim Area p1 p2 pn Eigen analysis of matrix C r provides shape information about the rim disk.Let us consider the vector λ = [λ 1 , λ 2 , λ 3 ] representing the eigenvalues sorted in decreasing order: λ 1 ≥ λ 2 ≥ λ 3 and [v 1 , v 2 , v 3 ] are the corresponding eigenvectors.The most dominant eigenvalues λ 1 and λ 2 represent Major and Minor axis of the rim disk respectively.The least dominant eigen direction will be the disk normal so n p = v 3 .Now, Rim Disk Area can be computed as:

(a) Top view
where n z = [0, 0, 1] and A 2D : R 2 → R is the area of the projection of A 3D on XY-plane.Average Rim Disk Diameter (d r ): Let us consider p θ i and p θ i +π are two rim points corresponding to two parameterized opposite radial scans(Q(t, θ i ) and Q(t, θ i + π)) then average rim diameter is written as: Average Pit Depth (h p ): represents the depth of foveal pit and can be computed as the difference between the average rim height and central foveal thickness: Average Maximum Pit Slope (s m ): Average maximum pit slope measures the steepness of foveal pit.It is defined as the average of the maximal slope of each radial scan.The slope of a parameterized radial scan Q(r, θ i ) is defined as: where Q y I (t) and Q x I (t) represent parameterized x and y coordinates for the interior segment.To calculate the maxima of the above equation, the gradient, ∇ t s i = 0 is computed.This gives the value of t m i for the maximal slope (s m i ): , where P y 3 and P y 0 are the y coordinates of the end control points of Q I (t, θ i ).Similarly, P x 3 and P x 0 are the x coordinates of the end control points of Q I (t, θ i ).Average maximum pit slope of a volume is defined as: Average Slope Disk Diameter (d s ): Average slope disk diameter is computed similar to average rim diameter d r .Slope width is computed between two opposite parameterized radial scans: Q(t, θ i ) and Q(t, θ i + π) and the corresponding maximum slope points are p s θ i and p s θ i +π such that: where t m represents maximum slope point in parametric domain as shown in Eq. ( 24) and Q x (t m , θ i ), Q y (t m , θ i ) are the corresponding x and y coordinate.Average slope width of a volume is defined as: Slope Disk Area (A s ): Slope disk area can be computed similar to the rim disk area.Let us consider (p s 1 , p s 2 , • • • , p s n ) ∈ R 3 are maximum slope points corresponding to each radial scan and are computed using Eq. ( 24) and (26).Covariance of maximal slope points can be computed using Eq. ( 20) to obtain the normal vector, major and minor axes of the slope disk.Pit Flat Disk Area(A f ): Pit flat disk area measures the flatness of foveal pit around the center and is computed using a threshold value τ for each of radial scan.In each radial scan, a point p where retinal thickness is smaller than τ is computed.Then the corresponding segment from center (x m , y m ) to the computed point p f θ i is treated as flat.This can be done in following two steps: 1. First of all, the parametric value t f corresponding to threshold value τ is calculated: 2. Now, the x value corresponding the parameterized value t f is computed: Point p f θ i = (d s i , τ, θ i ) corresponds to the pit flat point for a single radial scan.Similarly, let us consider (p 3 as pit flat points.The area, major and minor axis of the pit flat disk can be computed similar to rim disk area.Average Pit Flat Disk Diameter (d f ) is computed using d f i from Eq. ( 29): Rim Volume (V r ): In general, the volume under a surface in Cartesian and Polar domains and can be computed using the following equation: After the discretization process of the whole volume into n radial directions, the above equation will become: The above equation is modified using the following substitution.The interior parameterized curve Q I (t, θ i ) is defined between center and rim points.So the parametric value t → {Q x I (t, θ i ), Q y I (t, θ i )} is t = 0 and t = 1 at center and rim points respectively.Similarly, dt.After the substitution, rim volume can be written as: In the above equation, R = P x 3 and at end point of the parametric curve t is equal to 1. R might have different values for different radial scans because of asymmetry of foveal pit.Inner Rim Volume (V I R ): Is defined as the volume of fovea within a fixed radius from center.The radius is given by the user.Let us consider t u to be the parametric value corresponding to the user input radius r u for each radial direction.Now, Eq. (31) can be modified as follows: scans of interior region of fovea.As mentioned in Eq. ( 5), α and β represent the distance between control points P 0 , P 1 and P 3 , P 2 respectively.The lowest ICC value obtained for β m indicates that the distance between the P 3 and P 2 is not consistent.Perhaps, the small inaccuracy in rim point detection is leading to the low repeatability of the β m .On the other hand, α depends on center or beginning control point(P 0 ) which is fix for a volume scan.Therefore, α m has a better repeatability compared to β m .

Model accuracy
For root-mean-square-error (RMSE) comparison, we have implemented two state-of-the-art methods for foveal shape analysis, the one proposed by Ding et al. [25], and the one described in Dubis et al. [20] and compared them with our method.For a better readability, we renamed [25] as M 1 and [20] as M 2 .In case of M 1 the 3D approach was implemented, as presented in the paper [25].For M 2 the 3D version was modeled by extending the 2D method presented by the authors on our radial re-sampled volume.In order to visualize the behaviour of M 1 vs. M 2 vs. CuBe, we present fitting results on the central B-scan for the same subject in Fig. 10.As shown in the figure, CuBe is able to capture the exact foveal shape with the lowest RMSE compared to methods M 1 and M 2 .Fig. 10.Comparison of the proposed method with two state-of-the-art methods [25] and [20].Blue curve is the raw input (ILM-RPE difference) at the central B-scan and black, red and green are the reconstructed curves with [20], [25] and our method (CuBe) respectively.The RMSE values show that the proposed method manages to reconstruct the shape with minimum artifacts, specifically the pit shape.
Methods M 1 and M2, both are using Gaussian-based basis function to model foveal shape.The method M 2 [20] uses DoG (difference of Gaussian) as basis function for each radial scan and derives only three parameters (slope, rim diameter and rim height) to characterize the shape of the fovea, which is a too small number to encounter all the variation of the fovea shape.The optimization procedure of the mentioned energy function in this method is complex because of the two different Gaussian basis and the energy function can have unstable minima quite often.The method M 1 [25] applies a Gaussian along with a quadratic and linear basis function to produce a 3D model of a volume scan.It needs eight parameters (A 0 , A 11 , • • • , A 22 ) to reconstruct a 3D model of a volume scan.However, these parameters do not represent any morphological information of the fovea, which makes their interpretation in relation to morphological characteristics rather difficult.Scheibe et al. introduced a stable and accurate modeling of the fovea using double derivative of an exponential function [22].This method [22] is applied to re-sampled radial scans from 3D macula cube scans.There are four parameters used to parametrize a radial scan.However, the parameter γ in the exponential basis function ( exp(−x γ ) [22]) makes the algorithm more complex and leads to difficulties in optimization and analytical analyses.Additionally, direct analytical derivation of 3D parameters, e.g volumes presents several computational problems.Recently, Liu et.al [24] has introduced a sloped piecemeal Gaussian model for characterizing foveal pit shape.This method uses a combination of linear and Gaussian basis function along with an additional parameter which encounters the pit bottom flatness.As this method is using piecewise basis function, the optimization process is not straight forward.However, this method only characterizes the foveal pit and not the complete fovea shape.
In Fig. 11, we show different RMSE values for 3D fitting on a volume scan.These RMSE values are computed between raw ILM-RPE segmented data (as mentioned in section 3.1) from Heidelberg Spectralis OCT and the corresponding fitted model.The comparison has been performed between the proposed method, CuBe and methods M 1 , M 2 .The RMSE values for our approach show overall lower values compared to other methods, as our method reconstructs the foveal pit more accurately as shown in Fig. 10.In the proposed method, the highest RMSE occurs in regions where the segmentation is mostly influenced by blood vessels.These blood vessels produce several "jumps" spatially close to each other as shown in Fig. 12 (bottom).These "jumps" have a strong influence also in the detection of maximum height points (rim points) which can lead to a higher RMSE because our parametrization scheme depends on these critical points.In such cases, peaks are representing the vessels and not the retinal tissue, which might induce an additional noise when investigating differences between healthy and pathological data [28].
An important aspect of the presented method is the ability to utilize characteristic fovea properties for the whole circular region (re-sampled radial region) because of the flexible and robust cubic Bezier parametrization scheme.This implies that we derive the discussed 3D parameters in a completely new fashion.This gives us the possibility to explore new morphological features of fovea in terms of Rim Disk, Slope Disk, and Pit Flat Disk by looking at the eigenvectors of the covariance matrix defined in Eq. (20).By computing pit flatness and its area from the 3D reconstruction, we are able to provide not only a better visual insight of this region but also a potentially new diagnostic parameter for further investigations into several diseases e.g.those characterized by nerve fibers and ganglion cells loss [4,7,8].Another beneficial aspect of the presented 3D approach is in the flexibility of defining the foveal radius.This can be computed from the model itself as a specific feature of each radial Bézier curve part, and as such it would take into account the asymmetric nature of this region, or as a variable radius, defined by the user.

Application in HC and patients with autoimmune neuroinflammatory disorders
Having all the parameters described in the previous section, a first analysis of the morphological aspects of the data set can be performed.Table 2 shows all the 3D parameters defined for the HC and patient groups.The measurements obtained with our method have similar values to the ones encountered in literature.For characteristics like Rim Disk Area, Slope Disk Area, Pit Flat Disk Area a comparison to existing literature was not possible since this is the first time that these parameters have been introduced.
We also investigated the capability of the derived parameters to differentiate between HC and patients.To this end GEE analyses was performed.Table 2 presents that several parameters show significant differences between HC and patients.These results open new possibilities for further investigation of fovea shape derived parameters in more specific clinical diseases.

Conclusion
In summary, we have developed a reliable, accurate and meaningful approach for fovea shape analysis, which is able to correctly model the profile of the foveal region and reconstruct its 3D shape.Our method has been shown to robustly encounter possible variations in foveal shape in HC but also in foveas that undergo considerable changes during the course of a neuroinflammatory disease.
The mathematical model created is simple and has the advantage of a straight forward derivation of parameters.The computed parameters are in direct relation to the geometry modeled, and therefore provide an intuitive way of interpretation for further medical analysis and clinical interpretation.A major advancement of the developed method is that it is possible to analyze the foveal shape in a clinical context, especially from the 3D perspective.
Several derived foveal shape parameters showed statistically, significant differences between HC and patients with neuroinflammatory diseases of the central nervous system and could potentially reveal more insights into the foveal morphology and the changes it undergoes especially when correlated with other clinical information.Further applications in ophthalmologic diseases like macular degeneration are worth to be investigated as well.
Table 2. Analysis of all the 3D parameters defined for the HC and patient group.The last column shows the GEE analysis between the two groups.Abbreviations: HC -healthy controls.SD -standard deviation, Min -minimum value, Max -maximum value, GEEgeneralized estimating equation models analysis accounting for the inter-eye/intra-subject dependencies, p -p value

Fig. 1 .
Fig. 1. Figure shows the structures of interest in the retina with layers ILM (Inner limiting membrane), GCL (Ganglion cell layer) and RPE (Retinal pigment epithelium layer).

Fig. 3 .
Fig. 3. (a) The central B-scan of a volume with invariant tangent at critical points.Red dots = critical points.(b) The interior and the exterior segment of the corresponding right half of the B-scan.The inner (P 1 , P 2 ) and end (P 0 , P 3 ) control points are shown in green and red, respectively.

Fig. 4 .
Fig. 4. The pipeline of the proposed algorithm.

3 .
Segmentation of the Radial Scans: Each of the radial scans c(r, θ i ) is segmented into interior c I and exterior c E regions at corresponding maximum (critical) point, as shown in Fig. 3(b).

Figure 5 (Fig. 5 . 5 .
Fig. 5. 3D shape reconstruction procedure.(a) Shows a volume scan with 61 B-scans and 768 A-scans and corresponding 24 radial directions (blue) using the bilinear interpolation.(b) Represents the 24 fitted radial scans (green) using the least square optimization.The red points show the critical points for each radial scan.

Fig. 6 .
Fig.6.2D shape reconstruction (a) Shows re-sampled radial scans from a volume scan divided in two parts: c(r, θ i ) and c(r, θ i + π) (b) Shows corresponding fitted radial scans using least square optimization.These are represented by Q(t, θ i ) and Q(t, θ i + π) and i = 1, • • • , n.A full radial scan can be consider as a combination of c(r, θ i ) and c(r, θ i + π) and is shown in different colors and corresponding parametrized curves are shown in the same color.

Fig. 8 .
Fig. 8. a: Top view of rim disk.b: Normal vector of rim disk plane and corresponding covariance eigenvalues and eigenvectors.

Fig. 9 .
Fig.9.A visual representation of the 3D parameters.The rim and the inner rim volume is defined as the volume covered by the corresponding disk area.The radius for the inner disk area and the volume is defined by user.

Fig. 12 .
Fig. 12.Two examples of 2D curve fitting results with the lowest (top) and the highest (bottom) RMSE values selected from the entire data set analyzed.Original data is shown in blue and fitted curve in green.RMSE values are from top to bottom, 1.0 µm, 7.1 µm.
The voxel dimensions in horizontal and, axial directions and distance between B-scans in this data set were approximately 11.69, 3.87, and 125µm respectively.A total of 187 OCT scans, 95 from healthy controls (HC) and 92 from patients with different autoimmune neuroinflammatory diseases, were selected from the NeuroCure Clinical Research Centers' imaging database.