**Purpose.**:
To use surfaces generated by two-dimensional penalized splines (2D P-splines) to characterize the shape of the macular ganglion cell plus inner plexiform layers (GCL+IPL) in a group of normal humans.

**Methods.**:
Macular images of the right eyes of 23 normal subjects ranging in age from 18 to 75 years were obtained with spectral-domain optical coherence tomography (SD-OCT). The thickness of GCL+IPL was determined by manual segmentation, areas with blood vessels were removed, and the resulting maps were fit by smooth surfaces in polar coordinates centered on the fovea.

**Results.**:
Smooth surfaces based on 2D P-splines could precisely represent GCL+IPL thickness data, with errors comparable to the axial resolution of the SD-OCT instrument. Metrics were developed for the size, shape, and slope of the edge of the foveal depression and size and shape of the surrounding macular ridge. The slope of the foveal edge was negatively correlated with foveal size (*r* = −0.60). The size of the macular ridge was positively correlated with foveal size (*r* = 0.75), with a slope near unity (0.90 ± 0.18). The centroids of the foveal edge and macular ridge clustered near the foveal center. The foveal edge and macular ridge were well fit by ellipses. The mean GCL+IPL thickness formed an elliptical annulus elongated by approximately 30% in the horizontal direction.

**Conclusions.**:
The methods developed here provide precise characterization of retinal layers for the study of glaucoma, foveal development, and other applications.

^{ 1 }Thus, the shape of the GCL in the normal macula becomes important background knowledge for the detection of glaucoma. The shape of the GCL is determined by the centrifugal displacement of the inner retina during formation of the fovea,

^{ 2–4 }so studying the shape of the GCL also provides insight into foveal development. Because the boundary between the GCL and the adjacent inner plexiform layer (IPL) is often indistinct, the thickness of the GCL and IPL combined (GCL+IPL) often serves as a surrogate for GCL thickness in clinical applications.

^{ 5–8 }

^{ 9 }A similar shape is evident in retinal thickness maps obtained by OCT in vivo.

^{ 1,10,11 }

^{ 12,13 }This method has been extended to use separate fits to multiple radial profiles,

^{ 14 }but is limited by the global properties of the particular model, which necessarily link central and peripheral features, and cannot reveal deviations from the chosen analytic form. Ideally, analysis of higher-dimensional data requires a compact mathematical representation with the ease of evaluation of an analytic function and enough flexibility to capture the essential structural variation of the features of interest.

^{ 15–17 }can be used to describe a thickness map as a smooth analytic surface while imposing minimal restrictions on its geometry. The P-spline method yields an array of coefficients that represent local properties of the original data and provides precise control of data smoothing in different dimensions as a means to reduce measurement noise. Additionally, the method can easily handle missing data and, when applied to a rectangular data array, the method is fast and computationally efficient.

^{ 16 }This study used 2D P-splines to quantify GCL+IPL thickness maps to determine the size and shape of the foveal depression, the size and shape of the surrounding macular ridge, and the relationships between them in a group of normal subjects.

^{ 18 }and who were normal according to its criteria were recruited to have extra OCT scans of their maculas. Data were from the right eyes of 23 subjects, 11 females and 12 males, with ages that ranged from 18 to 75 years and a median age of 48 years. Eyes were selected as described elsewhere

^{ 19 }to cover a wide range (>2:1) of foveal sizes. All subjects provided informed consent to participate in this study, which was approved by the Institutional Review Board of the University of Miami.

^{ 19 }

**Figure 1.**

**Figure 1.**

^{ 15 }Eilers et al.,

^{ 16 }and Currie et al.

^{ 17 }) for more detailed descriptions of the implementation.

^{ 20 }The fits in Cartesian coordinates used 25 order-four (cubic) B-splines in each direction with both directions penalized by λ = 0.6. The fits in polar coordinates used 16 order-five (quartic) B-splines in the radial direction with λ = 0.2 and 39 order-four (cubic) B-splines in the angular direction with λ = 3.0. Higher-order B-splines were used radially to provide more accurate radial derivatives.

^{ 21,22 }The number of basis functions was chosen to give a reasonable representation of the data. Computational time was quite fast (<0.5 seconds) and was not an important factor in the choice. After the original thickness data were fit, only the smooth surfaces were used in further analyses.

^{ 23 }Linear regressions were performed using a commercial spreadsheet application program (Excel; Microsoft Corp., Redmond, WA). Descriptive statistics are reported as mean ± 1SD.

**Figure 2.**

**Figure 2.**

**Figure 3.**

**Figure 3.**

^{ 9 }the GCL+IPL thickness in Figure 2 shows a deep fovea surrounded by a wide parafoveal (macular) ridge. As seen in the profiles of Figure 3, GCL+IPL thickness started near zero at the foveal center, rose steeply along the foveal edge (arrows) to the peak of the macular ridge, and then slowly fell. The rise along the vertical meridians (±90°; Figs. 3B, 3C) appeared steeper than the rise along the horizontal meridians (0°, 180°; Figs. 3A, 3D). The macular ridge was further from the center for the horizontal meridians, reflecting the oval shape of the contour lines seen in Figures 2B, 2C.

^{ 21,22 }To reduce the variation that might arise from using only a single point, the foveal edge was defined as all points with derivatives ≥75% of the maximum (marked by the red bar in Fig. 4A and the red segment in Fig. 4B), subject to the constraints that the points were ≥0.2 mm from the foveal center and that the maximum derivative was at a GCL+IPL thickness ≥30 μm. The slope of the foveal edge was determined by linear regression (Fig. 4B, dashed line).

**Figure 4.**

**Figure 4.**

**Figure 5.**

**Figure 5.**

*r*= −0.60,

*P*= 0.002).

**Figure 6.**

**Figure 6.**

**Figure 7.**

**Figure 7.**

**Figure 8.**

**Figure 8.**

*r*= 0.75,

*P*< 10

^{−4}). The regression line in Figure 9 has a slope near unity (0.90 ± 0.18).

**Figure 9.**

**Figure 9.**

**Figure 10.**

**Figure 10.**

^{ 24 }Seven examples selected to cover the range of macular and foveal sizes are shown in Figures 10A–G. The corresponding points are labeled in Figure 9. The SE of the fit for the foveal edge ellipses was 16.3 ± 6.4 μm, or approximately 3.7% of the average radius of the foveal edge. The SE of the fit for the macular ridge ellipses was 40 ± 8.7 μm, or approximately 3.6% of the average radius of the macular ridge.

*a*and

*b*are the lengths of the major and minor axes of the ellipse, respectively, and is the amount by which a circle must be stretched in the direction of the major axis to produce the ellipse. The tilt of an ellipse is the direction of its major axis. On average the macular ridge is more elliptically elongated (E

_{m}= 30.8 ± 9.3%) than the foveal edge (E

_{f}= 17.7 ± 9.3%).

^{ 23 }A map of the mean GCL+IPL thickness of the 23 right eyes studied is shown in Figure 11A.

**Figure 11.**

**Figure 11.**

^{ 9 }

**Table.**

**Table.**

Parameter | Foveal Edge | Macular Ridge |

Average radius, mm | 0.435 | 1.11 |

Centroid [x, y], μm | [–8.8, +3.5] | [–8.8, −4.0] |

Fitted ellipse | ||

SE of elliptical fit, μm | 8.4 | 22 |

Major semiaxis, mm | 0.468 | 1.27 |

Minor semiaxis, mm | 0.407 | 0.98 |

Tilt from nasal meridian | −2.5° | +3.7° |

Elliptical elongation | 14.9% | 29.0% |

^{ 9 }and OCT layer thickness

^{ 1,10,11 }and was exhibited by all 23 maculas studied. To characterize quantitatively the shape of the GCL+IPL, smooth surfaces generated by 2D P-splines in polar coordinates were used to define and measure two features: the foveal edge and the macular ridge. The locations of both features could be well fit by ellipses, which captured their tilt and elongation (Fig. 10H). The tilts of most of the foveal and all of the macular ellipses spread tightly around the horizontal meridian. The relative contributions to this spread of head position, ocular torsion, and anatomical variation are unknown and require further study. The greater elongation of the macular ridge was consistent with the shallower slope of the foveal edge along the horizontal meridian (Fig. 5A). The average radii of the foveal edge and macular ridge provided measures of their size that would be suitable for studies of demographic variation.

^{ 13,25,26 }

^{ 19 }and ganglion cell loss in glaucoma can be assessed with a thickness deviation map formed by subtracting the normal macular shape from a shifted version of a patient's map.

^{ 19 }

*r*= 0.009) and macular ridge (

*r*= 0.012), but the narrow distribution of axial lengths represented by the sample in this study precludes drawing any conclusions from this result. Future studies should include eyes selected to represent a wide range of axial lengths.

^{ 27 }

^{ 2–4,14,25 }The GCL+IPL represents a major fraction of the inner retina. The finding (Fig. 8) that the centroids of the foveal edge and macular ridge fall so close to one another and to the peak of foveal OS length suggests that under this hypothesis the forces involved in foveal development are rather accurately centered on a single point. In addition, the good correlation, with nearly unity slope, between the sizes of the foveal edge and macular ridge (Fig. 9) supports the idea that they move together during development. There are, however, individual differences in the relation, as seen in Figure 10. Including the GCL, IPL, and other retinal layers in studies of the relationship between the FAZ and foveal pit morphology

^{ 14 }may elucidate these differences.

*s*(

*x*,

*y*) can be formed by the summation of a set of appropriately chosen, localized 2D basis functions, that is, where

*B*(

_{ij}*x*,

*y*) are the basis functions and

*β*are coefficients that determine the shape of the surface. Each localized basis function is different from zero only in a relatively small portion of the domain where it is defined, but taken together the

_{ij}*β*cover the entire area. Given a set of basis functions

_{ij }*B*the task of fitting such a surface to an array of data can be accomplished by determining the

_{ij}*β*that minimize the error between the data and the surface. The basis functions used for 2D P-splines are the tensor products of B-splines (Fig. A1), smooth, one-dimensional (1D) curves formed from piecewise continuous segments of polynomials.

_{ij}^{ 21 }The B-splines used have uniformly spaced knots that span the range of the data grid. Forming each 2D basis function from two 1D functions greatly simplifies the fitting procedure.

**Figure A1.**

**Figure A1.**

**M**be an

*m*×

*n*matrix of data. The B-spline number and order (polynomial degree + 1) can be chosen independently along the two dimensions of the

*m*×

*n*grid. Each B-spline is evaluated at each point of the grid, so that the values for the B-splines in each direction form matrices, with one column for each spline. Following the notation of Eilers et al.,

^{16}let

**B̆**denote the

*m*×

*K*matrix of B-splines down the rows and

**B**the

*n*×

*L*matrix across the columns of

**M**, where

*K*and

*L*are the corresponding numbers of B-splines. The aim is then to fit

**M**with a surface written as where

**β**is a

*K*×

*L*matrix of basis function coefficients and

**B**′ denotes the transpose of

**B**. It should be noted that

**β**is typically much smaller than

**M**. Using equation 2, Eilers et al.

^{16}recently showed that, if a roughness penalty (from which the name P-spline originates) is used to control smoothing, the best fit is achieved when the matrix of coefficients

**β**satisfies the penalized regression equation The vec(·) function in equation 3 “flattens” a matrix into a column vector by stacking its columns, ⊙ is the element by element product of two matrices and

**W**is an

*m*×

*n*matrix of weights. The

*KL*×

*KL*matrix

**F**is a function of the spline matrices

**B̆**and

**B**and the weights

**W**. The

*KL*×

*KL*penalty matrix

**P**is described in the following text. This equation can be solved for vec(

**β**), a

*KL*× 1 column vector, which is then rearranged to yield the

*K*×

*L*coefficient matrix

**β**.

*m*×

*n*weight matrix

**W**allows one to use this approach even when some of the data are missing; the elements of

**W**that correspond to missing data are set to zero and the rest are set to one. The only constraint on

**W**is that the remaining data must be distributed over the grid in a way that provides adequate “support” for each estimated coefficient. It should be noted that the resulting fit is defined over the full grid and automatically extends the data over the regions where they are missing.

^{ 21 }With P-splines the roughness penalty is applied not to the curve itself, but to the coefficients that define the curve.

^{ 15 }The penalty matrix

**P**, an operator that calculates differences of the coefficients of the fit,

^{ 16,17 }is given by

**P**=

*λ*

**I**

*⊗ $ D \u2032 d $*

_{K}**D**

*+*

_{d}*λ̆*

**D̆**

*⊗*

_{D̆}**I**

*, where, respectively,*

_{L}**I**

*and*

_{K}**I**

*are identity matrices of sizes*

_{L}*K*and

*L*,

**D̆**

_{ D̆ }and

**D**

*are matrices formed by taking the*

_{d}*D̆-*and

*d*-th order differences along the rows of

**I**

*and*

_{K}**I**

*and*

_{L}, λ̆*λ*are penalty weights that can vary from 0 to ∞, and ⊗ denotes the Kronecker (tensor) product of two matrices. Note that the difference orders (

*D̆*,

*d*) and penalty weights (

*λ̆*,

*λ*) can be different for the two directions. Larger penalty weights correspond to more smoothing.

**x**and

**y**; for points on a grid

**x**and

**y**give the coordinates along the two edges of the grid; for points not on a grid

**x**and

**y**together form a list of coordinate pairs. The vectors

**x**and

**y**are used to calculate the B-spline matrices

**B**and

**B̆**, respectively, from the B-spline specifications. For points on a grid the surface is evaluated directly from equation 2 as

**S**(

*x*,

*y*) =

**B̆βB**′. For a list of points the surface is calculated from

**S**(

*x*,

*y*) =

**C**vec(

**β**), where

**C**= (

**B**⊗ $ e \u2032 K $) ⊙ ( $ e \u2032 L $ ⊗

**B̆)**is a matrix containing only the required tensor products of the B-splines and

**e**

*and*

_{K}**e**

*are column vectors of ones of lengths*

_{L}*K*and

*L*, respectively.

*Invest Ophthalmol Vis Sci*. 2011;52:2425–2436. [CrossRef] [PubMed]

*Vis Neurosci*. 2004;21:53–62. [PubMed]

*Vis Neurosci*. 2004;21:775–790. [PubMed]

*Vis Neurosci*. 2005;22:171–185. [CrossRef] [PubMed]

*IEEE Trans Med Imaging*. 2009;28:1436–1447. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2011;52:1412–1421. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2011;52:8323–8329. [CrossRef] [PubMed]

*Trans Vis Sci Tech*. 2012;1:3. [CrossRef]

*J Comp Neurol*. 1990;300:5–25. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2005;46:2012–2017. [CrossRef] [PubMed]

*Ophthalmology*. 2005;112:1734–1746. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2009;93:1223–1227. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2011;52:624–634. [CrossRef]

*Invest Ophthalmol Vis Sci*. 2012;53:1628–1636. [CrossRef] [PubMed]

*Stat Sci*. 1996;11:89–121. [CrossRef]

*Comput Statist Data Anal*. 2006;50:61–76. [CrossRef]

*J R Statist Soc B*. 2006;68:259–280. [CrossRef]

*Invest Ophthalmol Vis Sci*. 2011;52:7872–7879. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2012;53:3653–3661. [CrossRef] [PubMed]

*Functional Data Analysis*. 2nd ed. New York: Springer-Verlag; 2006:37–79.

*Functional Data Analysis with R and MATLAB*. New York: Springer-Verlag; 2009:29–44.

*Functional Data Analysis website*. Available at: http://ego.psych.mcgill.ca/misc/fda/. Downloads available at http://www.psych.mcgill.ca/misc/fda/downloads/FDAfuns/Matlab/. Accessed November 17, 2011.

*WCSG*. 1998. Available at: http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.1.7559. Accessed October 16, 2012.

*Invest Ophthalmol Vis Sci*. 2011;52:5105–5110. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2011;52:8769–8779. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2012;96:57–61. [CrossRef] [PubMed]