The Shape of the Ganglion Cell plus Inner Plexiform Layers of the Normal Human Macula

Abstract

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.

Introduction

The retinal ganglion cell (GC), the output neuron of the retina, sends its axon to the brain via the optic nerve. Disease, in particular glaucoma, can damage this axon, resulting in the death and disappearance of the GC. This structural loss can be especially apparent in the macula, where a significant fraction of GCs are concentrated to serve foveal vision and where the ganglion cell layer (GCL) is many cell bodies thick. Detection of GC loss for glaucoma diagnosis and management is potentially an important clinical application of optical coherence tomography (OCT). ¹ 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  

The overall shape of the GCL can be inferred qualitatively from maps of GC density obtained from flat-mounted human retinas, which show an elliptical mound elongated in the nasal–temporal direction with a steep-walled central depression surrounded by a ridge of high GC density that slopes down gradually toward the periphery. ⁹ A similar shape is evident in retinal thickness maps obtained by OCT in vivo. ^1,10,11  

Quantitative measurement of retinal features is greatly facilitated if their geometry is approximated by a smooth parametric surface; a complex set of data is replaced by a simpler system from which the measurements are extracted. For example, the macula has been modeled as a difference of Gaussians function to define measures of foveal size. ^12,13 This method has been extended to use separate fits to multiple radial profiles, ¹⁴ 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. 

The recently developed technique of two-dimensional penalized splines (2D P-splines) ^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. ¹⁶ 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. 

Methods

Subjects

Subjects who were participants in another study ¹⁸ 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 ¹⁹ 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. 

Macular Imaging and Layer Segmentation

Macular OCT imaging was performed with a commercial retinal imaging instrument (2D and 3D imaging) (Cirrus HD-OCT; Carl Zeiss Meditec, Inc., Dublin, CA). The 3D images contained 512 × 128 × 1024 voxels that sampled a nominal 6 × 6 × 2 mm region (horizontal × vertical × depth) centered on the macula. All images had good signal-to-noise ratio (signal strength >5) and those that exhibited eye movements during acquisition were excluded. Because the axial lengths of the study eyes were rather tightly distributed (23.8 ± 1.3 mm), no adjustment for axial length was used. The retinal layers were outlined by a combination of automatic and manual segmentations (Fig. 1). Automatic algorithms located the inner limiting membrane (ILM), the junction between the inner and outer segments (IS/OS), and the retinal pigment epithelium (RPE). The retinal nerve fiber layer (NFL) and IPL outer boundaries were manually segmented on selected B-scans using an interactive pen display (Cintiq 12WX; Wacom Technology Corp., Vancouver, WA). Where the NFL and IPL vanished in the fovea, their boundaries were made to coincide with the ILM. The location of the foveal center was selected at the point of maximum outer segment length (circle on IS/OS in Fig. 1). Retinal blood vessels were removed from analysis with a binary mask generated from a high-contrast OCT fundus image of the blood vessel shadows as described elsewhere. ¹⁹  

Figure 1. 

 

Layer segmentations (red lines) of a 6-mm-long horizontal B-scan through the fovea. ILM, inner limiting membrane; NFL, outer margin of the retinal nerve fiber layer; IPL, outer margin of the inner plexiform layer; IS/OS, inner edge of the junction between photoreceptor inner and outer segments; RPE, retinal pigment epithelium. The yellow circle indicates the foveal center.

Fitting GCL+IPL Thickness Maps

Because the method has not been applied previously to data from ophthalmic OCT, the Appendix briefly summarizes the formalism of 2D P-splines as used here. The interested reader should consult the references (Eilers and Marx, ¹⁵ Eilers et al., ¹⁶ and Currie et al. ¹⁷ ) for more detailed descriptions of the implementation. 

Surfaces can be fit to thickness maps in either Cartesian or polar coordinates. For Cartesian coordinates the surfaces were fit to data in the original 512 × 128 scan grid. For polar coordinates, data were converted by nearest-neighbor interpolation to a dense polar grid with its origin at the foveal center (point of maximum outer segment length; see Fig. 1) and the discontinuity at ±180° oriented horizontally (i.e., along the fast B-scan direction) on the temporal side. Barring large tilts of the subject's head during acquisition, the discontinuity will approximately follow the temporal meridian and temporal raphe. ²⁰ 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. 

All data fitting and most analyses were carried out using custom programs implemented with commercial technical computing software (Matlab; The MathWorks, Natick, MA). Matrices of B-splines and their derivatives were then calculated (Matlab function bsplineM.m, available on the Functional Data Analysis website). ²³ Linear regressions were performed using a commercial spreadsheet application program (Excel; Microsoft Corp., Redmond, WA). Descriptive statistics are reported as mean ± 1SD. 

Results

P-Spline Fits to a GCL+IPL Thickness Map

Figure 2 shows an example of data from one eye and the smooth surfaces that result from fitting with 2D P-splines. The GCL+IPL thickness in Figure 2A is the axial distance between the manually segmented boundaries labeled NFL and IPL in Figure 1. To reduce labor, only every other B-scan was segmented except near the fovea, resulting in the horizontal striped appearance. The mask formed from blood vessel shadows superimposed an additional branching pattern that contained no data. Figure 2B depicts the smooth surface that results from fitting the data with an array of 625 basis functions (25 × 25) in Cartesian coordinates. Figure 2C shows the same data converted by nearest-neighbor interpolation to a dense polar grid centered on the fovea and fit by an array of 624 basis functions (16 radial × 39 angular). Radial coordinates extended from 0.1 to 3.0 mm; angular coordinates extended from −180° to +180°, with 0° oriented nasally. 

Figure 2. 

 

Three representations of GCL+IPL thickness for a right eye. The yellow cross in each image marks the point of maximum OS length, which is assumed to be the foveal center. (A) Original thickness data from a 512 × 128 macular scan covering a 6 × 6 mm retinal area. Black areas contain no data. (B) A smooth surface fit to the data in Cartesian coordinates. The contour interval is 10 μm, with every other contour labeled. (C) A smooth surface fit to the data in polar coordinates centered on the fovea. Horizontal and vertical polar meridians are labeled at the image edge. The nasal and temporal edges of the image are marked N and T, respectively. Blue shading indicates areas of data that are not represented.

The macular geometry in normal subjects is expected to have a degree of symmetry around the fovea that favors the use of polar coordinates; structural features can be defined and measured along the radial dimension alone. The fits to the actual data in Cartesian coordinates then can serve to validate the fits to interpolated data in polar coordinates. The Cartesian and polar representations of Figures 2B, 2C were very close over their shared area (mean difference = 0.014 ± 1.2 μm). For the 23 eyes the average difference between Cartesian and polar representations was 0.03 μm, with an average SD of 1.4 μm. Polar coordinates were used for all analyses presented here. 

The relationship between the fitted P-spline surface and the data can be appreciated in Figure 3, which shows four profiles made along radial meridians extending from the foveal center in the eye of Figure 2. The data (black dots) were noisy, but the overall variation in GCL+IPL thickness is captured well by the fitted surface (red lines). For this eye, the standard error (SE) of the fit in Cartesian coordinates was 5.1 μm and in polar coordinates was 5.2 μm. For the 23 eyes studied, the SE of the fit in polar coordinates was 5.25 ± 0.38 μm. Thus, the scatter of the data around the smooth surface had a magnitude similar to the axial resolution of the OCT instrument. It should be noted that, with the penalties used, spline fitting did not smear the steeply rising edge of the fovea (Fig. 3, arrows). 

Figure 3. 

 

Profiles of GCL+IPL thickness along four radial meridians from the foveal center, two horizontal (A, D) and two vertical (B, C). The red lines show radial cuts through the smooth polar representation in Figure 2C and the black dots show the underlying data in Figure 2A. The arrows indicate the steeply rising foveal edge.

As expected from the known anatomy of the GCL, ⁹ 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. 

The slope and location of the foveal edge and the location of the macular ridge were defined for the local geometry of each radial profile and then global properties of these features were derived using a large number of equally spaced meridians. 

Foveal Edge

The slope of the rising foveal edge was determined from the region on the radial profiles where GCL+IPL thickness changed most rapidly. To find this region, radial derivatives of the smooth surface were evaluated from derivatives of the radial splines. ^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. 

 

Measurement of the rising edge of the fovea and the macular ridge on one meridian. (A) The derivative of the smooth profile in Figure 3D with the region defining the foveal edge marked by a red bar. (B) The profile from Figure 3D with the foveal edge marked in red. The slope of the foveal edge was defined by a linear fit to points on the edge (dashed line). The radial location of the 50-μm contour is marked with a dotted line and arrow. The parafoveal (macular) ridge is indicated by the horizontal black bar and its radial location by the black dot.

The slope and location of the foveal edge varied with meridian. Figure 5 shows the mean ± SD of three measures of the foveal edge as a function of meridian. The average slope was shallower horizontally (0°, ±180°) than elsewhere (Fig. 5A). The midpoint of the foveal edge might be used to characterize its location; the average GCL+IPL thickness at the midpoint (Fig. 5B) stayed within a fairly small interval (±5 μm) and was almost constant except near the temporal meridian (±180°). This suggested that a simple planar curve (in particular a thickness contour line) also might provide a useful definition of the location. The contour lines surrounding the foveal center (as in Figs. 2B, 2C) were approximately uniformly spaced, with the central 50-μm contour falling somewhat above the average thickness of the midpoint of the foveal edge (Fig. 5B, dashed line), so its radius was chosen to locate the foveal edge on each profile (Fig. 4B, dotted line with arrow). Figure 5C shows that the average radius of the foveal edge was largest temporally and smallest within superior and inferior sectors. Overall foveal size for an eye was defined as the mean radius of the 50-μm contour. 

Figure 5. 

 

Measures of the foveal edge as a function of meridian for 23 normal right eyes. The solid lines are population means; the dotted lines show ±1 SD. (A) Slope. (B) Thickness at the midpoint. The dashed line is at 50 μm. (C) Radius of the central 50-μm contour line.

Across subjects there was a negative correlation between the mean slope of the foveal edge and the foveal size (Fig. 6), that is, the edges of larger foveas tended to have shallower slopes (r = −0.60, P = 0.002). 

Figure 6. 

 

The slope of the foveal edge was negatively correlated with foveal size. Each point is one eye. The dashed line was fit by linear regression.

Macular Ridge

The wide macular ridge that surrounds the foveal depression was defined for each radial profile as the set of points that exceeded 90% of the profile peak (Fig. 4B). The macular ridge for the example eye is shown in Figure 7 as a translucent blue band. The location of the macular ridge (Fig. 4B, black dot) was defined as the thickness-weighted average of the radial positions of the points on the ridge and is shown in Figure 7 as the dotted red line, with one dot for each of 72 meridians spaced 5° apart. The location of the foveal edge, corresponding to the 50-μm contour line, is also shown (solid red line). 

Figure 7. 

 

Central 4 × 4 mm of the GCL+IPL thickness surface in Figure 2C with foveal and macular features superimposed. Solid red curve: location of the foveal edge. Pale blue band: the region that exceeds 90% of the peak for each radial profile. Dotted red curve: location of the macular ridge. Yellow cross: foveal center. Yellow circle: centroid of the foveal edge. Yellow dot: centroid of the macular ridge.

Figure 7 also shows the centroids of the foveal edge and macular ridge, which both lie close to the foveal center. The centroids of all 23 eyes clustered near the foveal center, as seen plotted on the central 100-μm square in Figure 8. The average distances between the centroids and the foveal center were 17 ± 7 μm for the foveal edge and 19 ± 12 μm for the macular ridge. The average distance between the two centroids was 20 ± 9 μm. 

Figure 8. 

 

The centroids of the foveal edge (blue) and macular ridge (red) relative to the foveal center (black cross).

The average radius across all meridians of the foveal edge and macular ridge, as defined above, provided convenient measures for the sizes of the fovea and macula, which for the 23 eyes were 0.44 ± 0.08 and 1.11 ± 0.10 mm, respectively. These measures were used to relate the two structures. As shown in Figure 9, the sizes of the foveal edge and macular ridge were correlated (r = 0.75, P < 10⁻⁴). The regression line in Figure 9 has a slope near unity (0.90 ± 0.18). 

Figure 9. 

 

The sizes of the foveal edge and macular ridge were positively correlated. Each symbol represents one eye and plots the averages over all meridians of the foveal and macular radii. The open circles labeled a–g correspond to the eyes in (A–G) of Figure 10. The × denotes the eye shown in Figures 2 and 7. The dashed line was fit by linear regression.

Figure 10. 

 

Ellipses fit to the shapes of the foveal edge and macular ridge. In (A–G) the scale is given by the 1 mm calibration bar in (G) and E_f, E_m are the elongation factors for the foveal and macular ellipses, respectively. The small dots are values calculated from the GCL+IPL surfaces and the red lines are fitted ellipses. For clarity only half of the calculated points are shown for the foveal edge. The foveal center is denoted by a cross. (A) Largest macula; (B) medium fovea, large macula; (C) medium fovea, medium macula; (D) medium fovea, small macula; (E) smallest macula; (F) smallest fovea; (G) largest fovea, which also had an unusual vertical orientation. (H) Elongation factors and tilt direction for all 23 foveal edges and macular ridges. For each ellipse a line extends in both directions from the origin in the direction of tilt. The length of this line is equal to twice the elongation factor. The arrow indicates the nearly vertically oriented ellipse fitted to the foveal edge of the eye in (G).

The oblong shapes of the fovea and macula could be well described by ellipses. The locations of the foveal edge and macular ridge were determined on 72 radial profiles and ellipses were fit to the resulting curves. ²⁴ 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. 

The properties of the foveal and macular ellipses for all 23 eyes are depicted in Figure 10H as lines that show the elongation factor and tilt of each ellipse. The elongation factor of an ellipse is given by $\text{E}=(a/b-1)\times 100\%,$ where 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%). 

For most foveas and all maculas studied, the ellipse was elongated approximately horizontally (Fig. 10H). The obvious foveal exception is shown in Figure 10G; although the macular ellipse of this eye was approximately horizontal, the foveal ellipse was nearly vertical, as indicated by the arrow in Figure 10H. The other foveal edge with a nearly vertical tilt angle had a very small elongation factor (4.3%), so that it differed only slightly from a circle (for which the tilt angle is undefined). 

Average GCL+IPL Thickness Map

The overall similarities of the macular GCL+IPL thickness among eyes was examined by forming a map of their mean. This was easily produced from the coefficients of the fitted surfaces in polar coordinates, because the maps for individual eyes were already centered on the fovea and the coefficients for the mean of a set of surfaces is simply the means of their coefficients. ²³ A map of the mean GCL+IPL thickness of the 23 right eyes studied is shown in Figure 11A. 

Figure 11. 

 

The mean GCL+IPL thickness of 23 normal right eyes. (A) Mean map created by averaging the smooth fits in polar coordinates. The contour interval is 10 μm, with every other contour labeled. The dashed square outlines the enlarged area in (B). (B) The central 4 × 4 mm of the mean map showing the locations of the foveal edge, macular ridge, and their centroids as in Figure 7. Black lines: best-fitting ellipses.

The mean thickness map had an approximately elliptical symmetry centered on the fovea, as shown by the contour lines in Figure 11A. As in individual eyes, with increasing radial distance from the fovea, the mean GCL+IPL thickness increased to a peak and then slowly decreased. The peak was approximately 10% thinner temporally than elsewhere, but otherwise the mean thickness map was remarkable for its uniformity. 

The locations of the foveal edge and macular ridge of the mean thickness map were both smooth curves, as seen in the expanded image of Figure 11B. The centroids of these curves fell very close to the common foveal center of the 23 eyes (Table). The shapes of the locations of the foveal edge and macular ridge in the mean thickness map were better fit by ellipses than for individual eyes, with SE of the fits approximately one-half the average individual values (Table). Both ellipses were nearly horizontal and, as for individual eyes, the macular ridge had approximately twice the elongation as that of the foveal edge. The elongation of the macular ridge found here with OCT was almost the same as that of the peak of the macular GC density found histologically (29% vs. 28%). ⁹  

Table. 

 

Foveal and Macular Parameters for Mean GCL+IPL Thickness

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%

Discussion

The shape of the GCL+IPL of the normal human macula is exemplified by the mean map in Figure 11; the GCL+IPL forms an elliptical mound in which a deep fovea with a steep edge is surrounded by a wide macular ridge that slopes gradually down toward the periphery. Qualitatively, this topography was expected from maps of GC density ⁹ 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  

These results have already shown potential clinical application; the overall symmetry in the shape of the GCL+IPL of the normal macula enables a new method for reducing variance for glaucoma diagnosis ¹⁹ 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. ¹⁹  

Essentially no correlation with axial length was found for the sizes of the foveal edge (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. ²⁷  

Although this study was not designed to address foveal development, it does provide some insight and suggests that the methods used may be useful to others. Theoretical modeling and anatomical data support the hypothesis that the forces of intraocular pressure and ocular growth act on the structural inhomogeneity of the foveal avascular zone (FAZ), to produce centrifugal displacement of the inner retina and centripetal displacement of the cones. ^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 ¹⁴ may elucidate these differences. 

To summarize, 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. Using this approach, 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. The size of the macular ridge was positively correlated with foveal size with a slope near unity. The centroids of the foveal edge and macular ridge clustered very close to 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. 

Acknowledgments

The authors thank Xiangrun Huang for providing the method used to form the blood vessel masks. 

Appendix: Method of 2D P-Splines

A smooth surface 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 β_ij 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 B_ij 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. ²¹ 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. 

 

Spatial relationship of the quantities used to fit a smooth surface to a data array. The red contours show one 2D basis function formed by the tensor product of two B-splines (red curves) along the two dimensions of the data array. The symbols along the rows and columns are used in the text. m, n: dimensions of the data array; K, L: dimensions of the coefficient array; B̆, B: B-spline matrices; D̆ _D̆ , D _d : difference operators of order D̆, d; λ̆, λ: penalty weights. In this illustration K = 11 and L = 13.

When working with data on a rectangular grid, the problem can be expressed in terms of matrices. Let 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.,¹⁶ 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.¹⁶ 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 β. 

The 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. 

The spatial variation of the data can be followed more accurately by a large number of small basis functions, but these may overfit the data and increase sensitivity to noise. Fewer, larger basis functions can provide a smoother fit, but may not accurately portray small details of a surface. This dilemma has been resolved by using smaller basis functions and then controlling the smoothness of a curve by penalizing its roughness. ²¹ With P-splines the roughness penalty is applied not to the curve itself, but to the coefficients that define the curve. ¹⁵ The penalty matrix P, an operator that calculates differences of the coefficients of the fit, ^16,17 is given by P = λ I _K ⊗ ${D^{\prime }}_d$ D _d + λ̆ Display Formula D̆ _D̆ ⊗ I _L , where, respectively, I _K and I _L are identity matrices of sizes K and L, D̆ _D̆ and D _d are matrices formed by taking the D̆- and d-th order differences along the rows of I _K and I _L, λ̆ and λ 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. 

An important advantage of 2D P-splines is the fact that all the information about the fit is contained in the matrix of coefficients and the specifications of the B-splines, in a manner independent from the grid on which the original data were defined. One can use equation 1 or its matrix equivalent equation 2 to evaluate the fitted surface on another set of points, for instance a different grid or a set of points describing a line or curve. The points are specified by two vectors 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^{\prime }}_K$) ⊙ ( ${e^{\prime }}_L$ ⊗ B̆) is a matrix containing only the required tensor products of the B-splines and e _K and e _L are column vectors of ones of lengths K and L, respectively. 

Footnotes

 Supported in part by National Eye Institute/National Institutes of Health Center Core Grant P30EY014801, Department of Defense Grant W81XWH-09‐1‐0675, and an unrestricted grant from Research to Prevent Blindness (New York).

Footnotes

 Disclosure: R.W. Knighton, Carl Zeiss Meditec, Inc. (C), P; G. Gregori, Carl Zeiss Meditec, Inc. (F), P