3D modeling to characterize lamina cribrosa surface and pore geometries using in vivo images from normal and glaucomatous eyes

En face adaptive optics scanning laser ophthalmoscope (AOSLO) images of the anterior lamina cribrosa surface (ALCS) represent a 2D projected view of a 3D laminar surface. Using spectral domain optical coherence tomography images acquired in living monkey eyes, a thin plate spline was used to model the ALCS in 3D. The 2D AOSLO images were registered and projected onto the 3D surface that was then tessellated into a triangular mesh to characterize differences in pore geometry between 2D and 3D images. Following 3D transformation of the anterior laminar surface in 11 normal eyes, mean pore area increased by 5.1 ± 2.0% with a minimal change in pore elongation (mean change = 0.0 ± 0.2%). These small changes were due to the relatively flat laminar surfaces inherent in normal eyes (mean radius of curvature = 3.0 ± 0.5 mm). The mean increase in pore area was larger following 3D transformation in 4 glaucomatous eyes (16.2 ± 6.0%) due to their more steeply curved laminar surfaces (mean radius of curvature = 1.3 ± 0.1 mm), while the change in pore elongation was comparable to that in normal eyes (−0.2 ± 2.0%). This 3D transformation and tessellation method can be used to better characterize and track 3D changes in laminar pore and surface geometries in glaucoma.


Introduction
Glaucoma is the leading cause of irreversible blindness worldwide [1][2][3] and results in the degeneration of retinal ganglion cell axons and the death of retinal ganglion cells. Several studies support the idea that the initial damage to these axons occurs at the lamina cribrosa, a 3-dimensional (3D) porous meshwork consisting primarily of flexible, collagenous beams located within the optic nerve head (ONH) [4][5][6][7][8]. Increases in intraocular pressure (IOP) in a monkey model of experimental glaucoma have been shown to result in a posterior bowing and thickening of the lamina cribrosa early in the disease process, potentially stretching the loadbearing laminar beams and altering laminar pore shape and geometry [9,10]. Changes in laminar pores could damage the encompassed retinal ganglion cell axons that pass through the laminar pores from the retina to the brain, potentially resulting in axonal degeneration and vision loss.
The lamina cribrosa has been studied in vivo in normal and glaucomatous eyes. Most reports have used spectral domain optical coherence tomography (SDOCT), a technique that provides cross-sectional views of the lamina with high axial resolution, to characterize changes in laminar surface depth with disease [8,11]. However, the poor lateral resolution inherent in commercial SDOCT systems (~15 μm) currently limits the suitability of this technique to monitor and quantify changes in laminar pore geometry. En face images of laminar pores have been acquired and analyzed in living eyes using fundus photography [12,13] and scanning laser ophthalmoscopy [14,15]. While these studies have improved understanding of laminar pore geometry in normal and diseased eyes, the size and total number of pores that could be resolved and the number of eyes in which laminar pores could be successfully imaged were limited due to the reduced lateral resolution imposed by uncorrected ocular aberrations. Adaptive optics has recently been used to correct for the eye's optical imperfections and yield high-resolution (~2-3 μm) in vivo images of laminar pores [16][17][18][19]. Studies using adaptive optics scanning laser ophthalmoscope (AOSLO) imaging [17][18][19] have dramatically improved the ability to quantify anterior laminar pore geometry, but have not accounted for the curvature of the anterior lamina cribrosa surface (ALCS), i.e., the fact that en face images acquired using an AOSLO represent a 2D projection of a 3D curved surface.
To obtain a more accurate representation of the native pore geometry in living eyes, we transformed 2D AOSLO images onto a 3D ALCS that was modeled using surface points marked on SDOCT images acquired through the optic nerve head. We computed the mean curvature of the 3D surface for each eye and compared pore parameters calculated from 2D and 3D transformed images taken in normal and glaucomatous eyes to assess the impact of 3D transformation on pore geometry.

Animal preparation
All animal care experimentation procedures were approved by the University of Houston's Institutional Animal Care and Use Committee and adhered to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research. The anterior lamina cribrosa was examined in 11 eyes from six normal rhesus monkeys (Macaca mulatta) between the ages of 1.6 -7 years. The anterior lamina cribrosa was also examined in 4 monkeys following the induction of unilateral experimental glaucoma (i.e., one glaucomatous eye with the fellow eye as a control) [20]. Monkeys were anesthetized with 20 to 25 mg/kg ketamine, 0.8 to 0.9 mg/kg xylazine, and 0.04 mg/kg atropine sulfate to minimize eye movements during in vivo imaging [21]. Monkey pupils were dilated using 2.5% phenylephrine and 1% tropicamide and a contact lens was placed on the eye to prevent corneal dehydration during imaging.

In vivo imaging
Wide field (20°) scanning laser ophthalmoscope (SLO) images of the ONH were acquired in all eyes using a Spectralis HRA + OCT (Heidelberg Engineering, Germany). En face SLO images were taken such that the anterior lamina cribrosa was at best focus. The same instrument was used to acquire 24-48 SDOCT cross sectional radial B-scans (20° in extent) centered on the ONH [ Fig. 1(a)]. The termination of the retinal pigment epithelium (RPE) / Bruch's membrane (BM) interface (defining Bruch's Membrane Opening [BMO]) and ALCS were manually marked on SDOCT B-scans [ Fig. 1(b)] using custom MATLAB code (The MathWorks, Inc., Natick, MA) to generate a point cloud in 3 dimensional space.
En face AOSLO reflectance images (840 nm) of the anterior lamina cribrosa were captured over a 1.5° field at a rate of 25 Hz during closed-loop adaptive optics correction [18,22]. The power of the 840 nm superluminescent diode (Superlum, Ireland) at the corneal plane was ~300 μW (which is more than 10 times below the maximum permissible exposure dictated by the ANSI standards for a 1.5° field size, 1.5 hour exposure and wavelength of 840 nm) [23,24]. Individual frames from the AOSLO videos were registered offline using a normalized cross-correlation technique to remove distortions caused by eye movements. Registered frames were averaged to produce images with high signal-to-noise ratios. A montage of the anterior laminar surface was created from registered AOSLO images [ Fig. 1(c)]. Laminar pore boundaries were manually marked using Photoshop (Adobe Systems, San Jose, CA) as previously described [18]. Pore area and elongation (ratio of the major to minor axis of an ellipse best-fit to the boundary of a given pore) were calculated in the en face montages. Because SLO, SDOCT, and AOSLO images were taken over fixed angular field sizes, the lateral distance covered by each scan on the ONH could vary depending on the eye's ocular biometry. We measured each eye's axial length, anterior chamber depth, and anterior corneal curvature using an ocular biometer (IOLMaster, Carl Zeiss Meditec, Dublin, CA) and incorporated these measurements into a model eye to convert the lateral field size in our images from angular (i.e., degrees) to physical (i.e., micrometers) units [18]. All SLO, SDOCT, and AOSLO images from a given eye were scaled to the same resolution (microns/pixel) and registered in each eye.

Three dimensional anterior lamina cribrosa surface modeling and characterization
A thin plate spline was fit to the marked ALCS point cloud to model the surface in 3 dimensions using the MATLAB curve fitting toolbox (The MathWorks, Inc., Natick, MA) [ Fig. 2(a)] [25,26]. The kernel used for the thin plate spline was 2 2 ( ) log( ) f r r r = (1) where f(r) is depth of the surface at a given radius, r, from the center of the kernel. The thin plate spline is an interpolator formed as the linear combination of multiple kernels, which minimizes the bending energy, E(f), defined as To estimate a smooth surface, a smoothing parameter (p) was added to the energy function for regularization [27]. The thin plate smoothing spline was estimated by minimizing the energy function, E s (f), defined as According to the MATLAB manual, the smoothing parameter, p, was estimated a priori. The smoothing parameter, p, was computed from the coefficient matrix of the linear equations obtained using a radial basis function applied to the (x,y) coordinates of the marked points: In our case, we found a consistently low value for the smoothing parameter (mean value of p = 1.6 x 10 −6 ± 6.1 x 10 −7 ) across all eyes. The standard deviation about the mean value was very small, indicating a very similar smoothing parameter across eyes.
Following this surface estimation, the mean anterior lamina cribrosa surface depth (ALCSD) was calculated as the mean perpendicular distance from a plane best-fit to the BMO point cloud to the thin plate spline surface (fit to the marked ALCS points) for locations within the BMO ellipse [ Fig. 2 The 3D anterior lamina cribrosa surface was also characterized by computing the mean radius of curvature of the surface for points located within the projection of the BMO ellipse onto the ALCS. In two dimensions, the radius of curvature (ρ) at a given point on the curve may be calculated as the ratio of the change in arc length (ds) to the change in the angle of a line tangent to the curve (dθ) between 2 points adjacent to the given point [Figs. 2(c), 2(d)] [28,29].
The curvature (κ) at the given point is the reciprocal of the radius of curvature: The mean curvature of the 3D laminar surface was calculated as the average of the values of the mean curvature calculated at points on the fitted surface (~1 μm separation) located within the projection of the BMO ellipse. The mean curvature at a given point on a 3D surface, H, can be calculated as the mean of the principal curvatures at that point. The principal curvatures (κ1, κ2) are the maximum and minimum curvatures at the given point: The mean radius of curvature for the 3D anterior laminar surface was calculated as the reciprocal of the mean 3D curvature. The AOSLO montage was then registered with the SLO and SDOCT images and projected onto the thin plate spline representing the ALCS [ Fig. 3(a)]. The 3D surface was tessellated into a mesh consisting of triangular elements with regularly sampled vertices in the x-y (lateral) plane spaced between 1.5 and 2 µm [ Fig. 3(b)]. The 3D area of a given pore was calculated by integrating the area of the tessellated triangles located within the pore boundary. If a triangle was not completely located inside the pore boundary, its contributing area was computed as the product between the ratio of the number of triangle vertices contained within the pore boundary and the triangle area. To compute the 3D elongation of a given pore, an ellipse was first fit in 3 dimensions to points located on the pore boundary of the 3D surface. The ratio of the major to minor axis of this best-fit ellipse was computed as the 3D elongation of the pore.

Validation of thin plate splines
Arbitrary known surfaces were created and used to validate and compare the goodness-of-fit of the thin plate spline with other fitting functions (e.g., bi-cubic spline). As shown in Fig. 4, thin plate and bi-cubic splines were fit to the same regularly sampled subset of points (25-100% of the total points) used to generate a known surface. The residual root-mean-square (RMS) errors of the thin plate and bi-cubic spline fits with respect to the known surfaces were calculated from the non-fitted points (or from all points if 100% of the points were used for the fit). For example, Zernike spherical aberration [30] was constructed as one of the known surfaces [ Fig. 4(a)]. Thin plate and bi-cubic splines were fit to 45% of the total known surface points that were regularly sampled. The RMS errors of the thin plate and bi-cubic splines (calculated from the residual 55% of the total known surface points) were 3.2% and 2.6% of the original peak-to-valley surface height of the known surface (3.35 μm), respectively [Figs. 4(b), 4(c)]. Similarly, we found the RMS errors for the thin plate and bi-cubic spline fits to be comparable when the data were regularly sampled for other known surfaces. However, because the ALCS points that were manually marked from the SDOCT images were not Fig. 4. (a) Zernike spherical aberration was created as a known surface for fitting validation. Surfaces were created by fitting a (b) thin plate spline and (c) bi-cubic spline to the same regularly sampled subset of points (45% of the total points) from the known surface. The RMS errors calculated from the non-fitted points were 3.2% and 2.6% of the total peak-to-valley height of the known surface (3.35 μm) for the thin plate spline and bi-cubic spline, respectively. regularly sampled (which is a requirement for a bi-cubic spline fit), we chose to use a thin plate spline as our fitting function (as further elaborated in the Discussion).
To validate the area measurements calculated using our tessellation method, the areas of surfaces with known parameters were compared with the areas generated after tessellating the known surface into a mesh of triangular planes. Figure 5 shows the area validation results obtained when using a hemisphere as a known surface. The difference between the theoretical area of the hemisphere and the area after tessellation was 0.6%. Similarly, the differences between the theoretical areas and the post-tessellated areas were less than 1% in all surfaces tested. These small differences were due to the accumulation of errors contributed by areas of the triangles not completely contained within the border (or edge) of the surface.

Three dimensional transformation of the lamina cribrosa in normal and glaucomatous eyes
Thin plate splines accurately fit the manually marked ALCS points (Fig. 6). The mean RMS fit error of the thin plate spline across the 11 normal eyes was 9.8 ± 1.2 μm, a value that is only 2.5 times larger than the SDOCT instrument's axial resolution (corresponding to a 2-3 pixel error). In addition, as seen in Fig. 6 and Table 1, the anterior laminar surface was relatively flat in normal eyes as exhibited by the large mean radius of curvature of the surface (mean = 3.0 ± 0.5 mm across eyes), and had a small mean ALCSD (190 ± 13 μm).
Despite having a more steeply curved ALCS [Figs. 7(c), 7(d)], the mean RMS error of the thin plate spline fit across the 4 glaucomatous eyes was 10.8 ± 1.3 μm, a value comparable to that measured in normal eyes. However, the mean radius of curvature of the laminar surface in the glaucomatous eyes was 1.3 ± 0.1 mm, demonstrating a more steeply curved (or "cupped") anterior laminar surface. This observation was also reflected in the larger mean ALCSD value (392 ± 31.4 μm) compared to normal eyes. The mean 2D pore area and elongation across all experimental glaucoma eyes were 1049 ± 389 μm 2 and 1.83 ± 0.64, respectively (Table 1). Following 3D transformation of the AOSLO images onto the thin plate spline surfaces, the mean pore area increased across the 4 glaucoma eyes by 16.2 ± 5.9% (160 ± 121 μm 2 ) due to the increased curvatures, while the mean change in pore elongation was −0.2 ± 2.0% (−0.005 ± 0.04). The mean 3D transformed pore area and elongation were 1209 ± 433 μm 2 and 1.84 ± 0.63, respectively (Table 1). While the anterior laminar surfaces in the glaucoma eyes were more steeply curved, the RMS fit errors for the thin plate splines were similar to those in normal eyes. The increased surface curvatures yielded a greater change in 3D pore area compared to that calculated in normal eyes.

Discussion
It is important to have a faithful representation of the anterior lamina cribrosa surface from the in vivo images as the ALCS serves as the basis for characterizing laminar surface and pore parameters in 3 dimensions. We considered parametric and non-parametric approaches to estimate the ALCS based on our manually marked data points. Parametric approaches, such as polynomial surfaces, were not suitable for estimating the ALCS because the shape of the marked data points did not conform to a single polynomial surface. For example, a lower order polynomial could provide a smooth surface fit to the marked data points in 3dimensional space with a potential consequence of yielding a high RMS fit error. Using a higher degree polynomial would yield a lower RMS fit error, but could result in undesirable ringing effects in the estimated surface (i.e., Runge's phenomenon) [31]. Hence, a nonparametric approach was found to be more suitable to describe the ALCS. Mutliple non-parametric approaches were considered for modeling the ALCS, including locally weighted regression, bi-cubic splines and thin-plate splines. Locally weighted regression surface estimation is a popular approach to estimate a smooth surface and model data points [32]. However, this smoothing surface approach requires a set of dense data points for regression within a local (or span) window that typically form a "point cloud" of data points which lie above/below and to the left/right of the actual surface. This type of "point cloud" is not generated in our data, as our data points are manually marked just along the top of the ALCS in each SDOCT B-scan. Therefore, a locally weighted regression approach was not suitable for our data.
We used a thin plate spline to model the anterior lamina cribrosa surface in 3D and quantified the goodness-of-fit of the thin plate spline method by calculating the residual RMS fitting error. This RMS error was found to be comparable to the residual RMS fitting error of a bi-cubic spline. Despite the fact that a bi-cubic spline tended to yield surfaces with slightly lower RMS fitting errors, we chose to fit all of our ALCS point clouds using thin plate splines because fitting with a bi-cubic spline requires that the data be regularly sampled, whereas the thin plate spline does not. The latter was more suitable for fitting a manually marked lamina cribrosa surface, as the generated point cloud was not regularly sampled [ Fig. 2(a)]. For example, it is often difficult to mark laminar surface points that reside underneath overlying blood vessels [see Fig. 1(b)]. Moreover, the residual RMS error was very low for the thin plate spline (Fig. 4). In addition, thin plate splines have been used previously to accurately represent the anterior and posterior laminar surfaces in histomorphometric reconstructions of the monkey lamina cribrosa [33] and have been shown to be an effective method for mapping and tracking non-rigid surface deformations longitudinally in the heart [34].
An accurate representation and characterization of the ALCS is necessary to better understand and track structural changes in the ALCS during disease progression. Several studies have shown an increase in mean ALCSD during the progression of glaucoma [8][9][10]. However, the use of this single parameter provides only a partial characterization of the structural changes in the laminar surface. While an increase in the mean ALCSD indicates an overall posterior shift in the anterior laminar surface (relative to the BMO reference plane), this value provides no characterization of the shape or curvature of the actual surface. For example, it is possible to have two differently curved anterior laminar surfaces with the same mean ALCSD; for instance, one surface could be uniformly depressed in all locations (i.e., a large radius of curvature) whereas the second surface could have a large degree of bowing or "cupping" (i.e., a small radius of curvature) about the same mean depth. Likewise, two anterior laminar surfaces with the same mean curvature could give rise to different mean ALCSD values if, for example, one surface is more posterior to the BMO reference plane than the other. Therefore, we believe that specifying the mean ALCSD and mean radius of curvature provides increased information on the structural changes that occur in the ALCS during disease progression and will facilitate better tracking of the laminar surface over time (Fig. 8). An increase in mean ALCSD accompanied by a decrease in the mean radius of curvature would indicate a more steeply curved laminar surface that is more posteriorly located with respect to the BMO reference plane, whereas an increased mean ALCSD accompanied by no change in the mean radius of curvature would indicate a more uniform posterior deformation of the entire surface. The 3D transformed pore parameters better represent the native state of lamina cribrosa surface pore geometry compared with values calculated from 2D AOSLO images. Table 1 shows that larger increases in the 3D transformed pore area were calculated in the experimental glaucoma eyes. These increases were due to the fact that the mean radius of curvature was lower (i.e., the curvature was greater) in glaucoma eyes relative to normal eyes ( Table 1). The small mean increases in laminar pore parameters calculated after 3D transformation in normal eyes were due to the fact that the anterior laminar surface was relatively flat with a high radius of curvature (see Fig. 6). However, the change in mean elongation of the laminar pores following 3D transformation was nearly identical between normal eyes and eyes with experimental glaucoma. The change in the elongation of a pore after 3D transformation depends on the spatial location and orientation of the pore. The 3D elongation of a transformed pore would increase relative to the 2D elongation of the same pore if the slope along the major axis of the ellipse best-fit to the transformed pore is greater than the slope along the minor axis of the same ellipse. Alternatively, the 3D elongation of a transformed pore would decrease from its 2D elongation value if the slope is greater along the minor axis of the ellipse best-fit to the pore. The fact that the mean changes in elongation following 3D transformation in normal and glaucoma eyes were 0.0 ± 0.2% and −0.2 ± 2.0%, respectively, of the original 2D pore elongation likely indicates the absence of a predominant orientation of the anterior laminar surface pores that were imaged.

Conclusion
As evidenced by the low RMS fitting errors, thin plate splines can accurately fit laminar surfaces with different slopes in normal and glaucomatous eyes. Changes in mean ALCSD and the mean radius of curvature of the anterior laminar surface can better describe structural changes that occur during disease. The laminar surfaces in normal eyes are fairly flat (i.e., exhibit large radii of curvature), yielding small changes in 3D pore geometry. As the curvature of the anterior laminar surface increases, which typically occurs during disease progression, larger changes can be noted in the 3D transformed pore geometry. The proposed 3D transformation and tessellation method can be employed with en face and axial imaging techniques, such as AOSLO and SDOCT imaging, to better understand 3D laminar geometry. This approach could be used in vivo to further examine global and regional differences in the 3D laminar geometry of normal eyes, 3D changes in laminar geometry that occur at different stages of the disease in experimental models of glaucoma and in human glaucoma patients, and 3D differences in laminar geometry between normal and glaucomatous eyes.