Automated segmentation of anterior lamina cribrosa surface: How the lamina cribrosa responds to intraocular pressure change in glaucoma eyes?

The lamina cribrosa is a sieve-like structure where retinal ganglion cell axons and central retinal vessels exit from the eye through the scleral canal. The lamina cribrosa has been known to play an essential role in the physiopathology of glaucoma and has been investigated as a potential location to identify early glaucomatous damage. Many researchers have been studying how the lamina cribrosa responds to intraocular pressure change, leading to axonal insults. Recently, 3D spectral domain optical coherence tomography (SD-OCT) Enhanced depth imaging (EDI), an optical imaging technique, has been proposed to improve OCT imaging of deeper retinal structures such as the choroid and the lamina cribrosa. However, the shadowing from vasculature and other reflective structures make the segmentation of the anterior lamina cribrosa surface difficult. In this paper, we present a new approach for the segmentation of the anterior lamina cribrosa surface. To deal with the complexity of the surface segmentation, we propose the use of a shape-constrained surface evolution method where the surface is refined iteratively using a non-local Markov random field based segmentation. The estimation of the model parameter is addressed using a Metropolis-Hastings algorithm. Our experiments showed a significant correlation between change in the intraocular pressure level and change in the position of the lamina cribrosa over time.


INTRODUCTION
Glaucoma, the second leading cause of blindness in the world [1], is a group of neurodegenerative eye diseases that leads to the loss of peripheral vision and, in an advanced state, irreversible blindness. Therefore, it is important to develop sophisticated tool for early glaucoma diagnosis in order to avoid permanent damage to the optic nerve head by slowing the progression of the disease. The primary damage, leading to a loss of vision, occurs to the retinal ganglion cell (RGC) axon. Many studies suggest that the lamina cribrosa (LC), a sieve-like structure that provides structural and functional support to the RGC axons, is the principal site of RGC axon insult and injury in glaucoma [2]. For many years, researchers used histological techniques to elucidate the details of LC structure, both in normal and glaucoma eyes. However, advances in optical coherence tomography (OCT) have provided a noninvasive optical imaging technique that has been used to visualize optic nerve head ONH structures. Recently, the Spectralis spectral domain SD-OCT Enhanced depth imaging (EDI) mode was proposed to image deeper retinal structures such as the choroid and the LC in order to improve glaucoma diagnosis and risk management [3]. The ability to visualize the LC in vivo allows us to better understand the disease. In particular, since the only treatment to slow glaucoma progression is by lowering intraocular pressure (IOP), regardless of the level of IOP, it is important to understand the effects of IOP on the LC.
Histologic studies of humans and animals as well as experimental animal models based on IOP elevation have shown that eyes with glaucoma often have deformities of the LC in particular posterior LC displacement [4]. Several studies propose the anterior lamina cribrosa surface depth (ALCSD) measurement to model such displacement [5]. The ALCSD parameter is defined as the average perpendicular distance from the anterior LC surface (i.e; the surface separating the prelaminar tissue and LC structure) relative to Bruch's membrane opening (BMO) (i.e; the termination of the Bruch's membrane (BM) layer) reference plane. Hence, ALCSD parameter requires the identification of the BMO and the segmentation of anterior LC surface.
A number of studies have tried to address the ONH layer segmentation problem. The 3-D graph search approach was used in [6] to segment ONH layer where the gradient magnitude of the dark-to-bright and bright-to-dark transitions was used as cost function. However, poor quality and noisy images, common in clinical management of glaucoma, may affect the layer segmentation. To overcome this problem, we proposed in [7] a new deconvolution-based approach for ONH layer segmentation to identify the BMO. However, the shadowing from vasculature and other reflective structures makes the deconvolution based approach no longer appropriate to segment the anterior LC surface. Due to these limitations, researchers rely on manual segmentation of the anterior LC surface which is a time-consuming process and impractical for routine clinical use.
In this paper, we propose a new Bayesian method to segment the anterior LC surface. Bayesian methods are relatively simple yet efficient and offer elegant tools to include priori knowledge in the inference model. Prior knowledge can be incorporated in the model such as information about the shape, the position and distribution of the LC voxels. In this work, we propose the use of the Markov Random Field (MRF) to exploit the statistical correlation of intensities among the neighborhood voxels. Recently, the MRF prior based methods have been used for medical image segmentation [8]. However, due to shadowing from vasculature and other reflective structures, the intensities of the voxels across the LC surface may have different appearance. To overcome this problem, we consider the use of a non-local framework for the MRF energy function definition which has been first proposed for image denoising tasks [9]. The main idea of the non-local approach is to exploit repetitive structures in the image which leads to a multi-model approach using only a single observation. The estimation of the model parameter is addressed using a Metropolis-Hastings algorithm [10]. Specifically, the sampling approach was inspired from the Biased and Filtered Point Sampling (BFPS) method [11] where authors proposed a perturbation-based sampling approach: to generate a new boundary of the object to segment, it is sufficient to perturb the previous boundary of the object. Therefore, the evolution of the segmentation is addressed iteratively by changing the surface of the LC structure according to the non-local MRF energy function. This approach differs from [11] in term of the segmentation strategy. While, in [11], authors only use the boundary information to segment the images, we combine in this method both the boundary information and the object information to improve segmentation of anterior lamina surface.
The paper is divided into two sections. In section 2, the proposed anterior LC surface segmentation is presented. Then, in section 3, we first compare the segmentation accuracy of the proposed scheme with the BFPS method [11] and then we study the effect of the change in IOP on change in LC position through the ALCSD parameter.

ANTERIOR LAMINA CIBROSA SURFACE SEGMENTATION
In this work, we utilize 3D Spectralis SD-OCT images (Heidelberg Engineering) radial scan images. Each 3D image consists of 48 enhanced depth imaging radial 2D B-scans centered on the optic nerve head. Each B-scan consists of several inter-retinal layers/structures (e.g; the BM layer, the Retinal Nerve Fiber layer (RNFL), LC structure). The aim of this work is to study the response of the LC to the IOP variation through the ALCSD parameter. To this end, we are only interested in the identification of the BMO points and the segmentation of the anterior LC surface. Figure 1 presents a B-scan image with the different parameters to estimate. The BMO points are identified using the method presented in [7]. The BMO points are then used to 1) estimate the BMO reference plane and 2) to estimate the area of interest (AOI) (c.f; figure 2) where the LC is localized in order to optimize the segmentation process (i.e; the LC is localized below the BMO reference plane). We denote by where Z is the normalization constant, Θ consists of the model hyperparameters and U (H j |Y k , Θ), j ∈ {0, 1, 2}, is the energy function of the MRF model. In this work, the isotropic Potts model-type prior is adopted, the energy function to minimize is: where p l (Y k (l)|Q(l) = H j ; Γ j ) is probability density function of the lth intensity Y k (l) belonging to the source neighboring system N i given the class H j , Γ j is the set of the distribution hyperparameters of the class H j , j ∈ {0, 1, 2}, δ the delta Kroneker function; and ζ is a positive parameter. As in [7], we assume that p l (Y k (l)|Q(l) = H 0,1 ; Γ 0,1 ) follows a Gamma distribution and p l (Y k (l)|Q(l) = H 2 ; Γ 2 ) follows a Gaussian distribution. Note that we opted for the 3D 8connectivity neighboring system. We propose a non-local approach to define the energy function. Because neighborhood around a voxel may be similar to the neighborhoods around other voxels in the same image, the non-local approach allows us to exploit this redundancy by extending the neighboring system. Therefore, the segmentation of a voxel will depend not only on its neighborhood but also on other neighborhoods called patches. The contribution of each new patch depends on the similarity between the main neighborhood N i and the other patches N i d where i d is the center of the patch. To measure the similarity of the patches, a weighted graph G that links together voxels over the image is estimated. The Euclidean distance has beem shown to be is sufficient to reliably measure the similarity [12]. Therefore, the similarity between two voxels i and i 1 , based on the pairwise Euclidean distance, is expressed as: The new energy function which integrates the non-local weights G associated with different patches around the voxel of interest i and located in an extended neighborhood N N i is expressed as: The first part of the energy function models the a priori we have on the AOI Y k (l) conditioned to each class as well as the information available in other patches (information fusion). The second part models the spatial dependency by the use of the second-order isotropic Potts model with parameter ζ. Hence, the definition of the energy function favors the generation of homogeneous areas reducing the impact of the speckle noise, which could affect the segmentation results. The hyperparameters ζ, η ∈ [0, 1] handle the importance of the energy terms. On one hand, ζ allows us to tune the importance of the spatial dependency constraint.On the other hand, η models the reliability of the patches. However, this energy function does not allow a control for the anterior LC surface shape. In [11], the authors proposed an elegant way to control for the shape of the segmented objects by sampling their boundaries using a Hastings-Metropolis sampler. To sample a new boundary, it is sufficient to deform the current boundary around a foot point by pushing it in a given direction using Gaussian filters. In this work, we opt for the same sampling approach (c.f; Figure 2). However, we used a digital topology approach to perturb the surface [13]. In order to ensure that each step satisfies detailed balance, the set of the predefined perturbations is designed such that the reverse to the previous step always exists. Hence, to sample the segmentation map Q, we first sample the surface S then we sample Q according to S and the observation Y k . The main Hastings-Metropolis sampler steps are described in Algorithm. 1. Moreover, in order to favor smooth surfaces, a curve length penalty param- Fig. 2. Description of the AOI and the disturbing-based approach for surface sampling on a B-scan (a radial cut of the 3D volume image) eter is added to the energy function: where β, ∈ [0, 1] handles the importance of the smoothness term. The segmentation map Q is estimated using the maxi- ) v) Until Convergence criterion is satisfied mum a posteriori MAP estimator.

RESULTS
The proposed anterior LC surface segmentation method called the non-local shape-constrained (NLSC) method was experimentally validated with clinical datasets. Eligible participants were recruited from the University of California, San Diego (UCSD) Diagnostic Innovations in Glaucoma Study (DIGS). In order to evaluate the anterior LC surface segmentation scheme, two experts manually marked the position of the anterior LC surface on 8 predefined vertical lines at each B-scan of 50 normal and 50 glaucomatous eyes. When the anterior LC surface is hard to see on one line, the expert does not mark the position (c.f; Figure 3 (a)). The failure rate (the number of marks of the anterior LC surface that are different from the reference) ranged from 0 to 2, by definition (0 when the mean difference < 3 pixels, 1 when the mean difference < 5 pixels and 2 when the mean difference > 5 pixels). Table 1 presents the different failure rates. We found that the NLSC method performs highly close results to both experts manual segmentation and perform better than the BFPS method.
To study the effect of the IOP on the LC structure (longitudinal study), we used 21 eyes from 21 glaucoma pa-tients followed longitudinally for 2.4 (1.6) years (median (interquartile range)). For each eye, we calculated the ALCSD parameter using the NLSC method. We used the R-squared coefficient to assess the correlation between the change in ALCSD and the change in IOP parameters over time. We found a significant correlation (R 2 =0.68) between the IOP variation and the ALCSD variation over time. Figure 4 (a) presents an example of a high correlated variation and Figure  4 (b) a moderate correlated variation of IOP and the ALCSD parameters. For the first case (Figure 4 (a)) the LC was deformed posteriorly or anteriorly according the IOP variation. This is not the case in Figure 4 (b) where the LC was more rigid to the IOP variation. More analyses is need to study this variability between patients which can constitute a big step forward in understanding the cause of the glaucoma disease and to develop sophisticated tool for early diagnosis.

CONCLUSION
In this paper, a new shape-constrained surface evolution method for the segmentation of the anterior LC surface has been developed. We particularly propose a non local approach to define the MRF energy function. The validation of the proposed approach, the first automated method to segment the anterior LC surface to our knowledge, on clinical dataset has shown its ability to perform similar results to the manual expert segmentation. Our experiments showed a significant correlation between the change in intraocular pressure and the change in position of the lamina cribrosa.