A Bayesian framework for glaucoma progression detection using Heidelberg Retina Tomograph

Glaucoma, the second leading cause of blindness
in the United States, is an ocular disease characterized by
structural changes of the optic nerve head (ONH) and changes
in visual function. Therefore, early detection is of high importance
to preserve remaining visual function. In this context,
the Heidelberg Retina Tomograph (HRT), a confocal scanning
laser tomograph, is widely used as a research tool as well as
a clinical diagnostic tool for imaging the optic nerve head to
detect glaucoma and monitor its progression.
In this paper, a glaucoma progression detection technique is
proposed using the HRT images. Contrary to the existing
methods that do not integrate the spatial pixel dependency in
the change detection map, we propose the use of the Markov
Random Field (MRF) to handle a such dependency. In order
to estimate the model parameters, a Monte Carlo Markov
Chain procedure is used. We then compared the diagnostic
performance of the proposed framework to existing methods
of glaucoma progression detection.


I. INTRUDUCTION
Glaucoma refers to a set of eye conditions of great clinical and etiological heterogeneity.It is characterized by the degeneration of optic nerve fibers and often accompanied by an elevated intraocular pressure (IOP).The loss of nerve fibers leads a decrease in the thickness of the retinal nerve fiber layer (RNFL), affects the appearance of the ONH and causes an irreversible damage to the retina.In the course of the disease, the neuroretinal rim gets thinner whereas the optic cup gets bigger [1].
Since the introduction of the ophthalmoscope by Helmholtz in 1851, ophthalmologists have been able to assess the ONH structure associated with glaucoma.Nevertheless, the qualitative clinical assessment of the ONH leads to a considerable inter-observer diagnostic variability [2].Therefore, development of quantitative measurements for glaucoma detection is important to make the qualitative clinical assessment more objective and reproducible.
For these reasons, sophisticated ocular imaging instruments are providing quantitative parameters of the ONH in glaucoma such as the Scanning Laser Polarimetry and the Optical Coherence Tomography.In particular, the Heidelberg Retina Tomograph (HRT; Heidelberg Engineering, Heidelberg, Germany), a confocal scanning laser tomography, has been commonly used for glaucoma diagnosis since its commercialization 20 years ago [3].
A limited number of studies have been published investigating glaucoma progression detection using HRT images by detecting changes between baseline reference images and follow-up images.In [4], authors assessed the glaucomatous changes using the Topographic Change Analysis (TCA).However, this method requires up to three additional confirmatory follow-up exams to detect changes (progression) [5].To overcome this requirement, [6] used the Proper Orthogonal Decomposition (POD) for glaucoma detection.However, although this method is successfully applied to HRT images, it does not exploit additional available knowledge, such as spatial dependency among neighboring pixels ( i.e, the status of a pixel will depend on the status of its neighborhood).The POD method indirectly utilizes the spatial relationship among pixels by controlling the familywise type I error rate.In [7], a retinal surface reflectance model and a homomorphic filter were used to detect progression from changes in the retinal reflectance over time.Similar to POD, dependency among spatial locations were only indirectly used by controlling familywise type I error rate.
In this paper, we propose a new strategy for glaucoma progression detection.We particularly show that glaucoma detection can be improved if spatial dependency of status of pixels are directly modeled and integrated within the change detection method.In addition, since for each follow-up exam, three follow-up scans are generated by the machine, each follow-up scan is considered as a separate 'information source'.The data-fusion (i.e.combining information from three scans per exam) and dependency in the status of pixels in a neighborhood are jointly modeled and addressed using the Markov Random Field (MRF) [8].Indeed, many studies have tackled the change detection problem by modeling the change detection map as a MRF [9], [10].The widely used procedure to estimate the different problem paramewww.ijacsa.thesai.orgters is the Expectation-Maximization EM algorithm [11].However, since we used the MRF model with the change detection map as a priori for the change detection map, the optimization step is intractable and therefore, we used a Monte Carlo Markov Chain MCMC technique [12] at each EM iteration [13].The principle of MCMC technique is to generate samples by drawing from the posterior densities [12].
The paper is divided into two sections.In section II, we present our new glaucoma progression detection scheme.In Section III, we present the results of our new algorithm using HRT follow-up exams of participants in the UCSD Diagnostic Innovations in Glaucoma Study.Then, we compare the diagnostic accuracy, robustness and the efficiency of our novel approach to three existing progression detection approaches, topographic change analysis (TCA) [4], proper orthogonal decomposition (POD) [6] and the reflectance based method [7].

A. Illumination correction
The HRT images can be affected by inhomogeneous background (illuminance) and effect may be differing among follow-up exams due to curvature of the retina and differences in the angle of imaging the eye between exams [7].Although, this problem is not due to glaucoma, it would have an influence on the subsequent statistical analysis and quantitative estimates of glaucoma progression.To this end, a reflectance-based normalization step is performed.Assuming that the optic nerve head is Lambertian, each HRT image I can be simplified and formulated as a product I = L x F where F is the reflectance and L is the illuminance.Because illumination varies slowly over the field of view, the illuminance L is assumed to be contained in the low frequency component of the image.Therefore, the reflectance F can be approximated by the high frequency component of the image.The reflectance component F describes the surface reflectivity of the retina whereas the illumination component L models the light sources.The reflectance image can then be used as an input to the change detection algorithm [7], [14], [15].
Several methods have been proposed to solve the problem of reflectance and illumuiance estimation including the homomorphic filtering [16], the scale retinex algorithm [17]- [19] and the isotropic diffusion approach [20].The reflectance based method of detecting progression uses a homomorphic filter to estimate retinal reflectance [7].In our new glaucoma progression detection scheme presented in this study, we used the scale retinex approach.The retinex approach has been successfully utilized in many applications, including medical radiography [21], underwater photography [22] and weather images enhancement [23].

B. Change detection
Let us consider the detection of changes in a pair of amplitude images.Change-detection can be formulated as a binary hypothesis testing problem by denoting the "change" and "no-change" hypotheses as H 1 and H 0 , respectively.We denote by I 0 and I 1 two images acquired over the same scene at times t 0 and t 1 , respectively (t 0 > t 1 ), and coregistered.A difference image R with N pixels was estimated as pixelwise difference between images I 0 (a baseline image) and I 1 (a follow-up image): R = abs(I 0 − I 1 ).HRT acquires three scans during each exam.Therefore, we obtain three image differences R(l, i) = {r(l, i)|l = 1, 2, 3, i = 1, 2, ..., N}.The change detection is handled through the introduction of change class assignments Q = {q(i)|i = 1, 2, ..., N}.Accordingly, the posterior probability distribution of change (i.e.q(i)) at each pixel location is expressed as: where Z is the normalization constant, Θ consists of the model hyperparameters and U (H j |R, Θ), j ∈ {1, 2}, is the energy function of the MRF model.The energy function is expressed as a linear combination of elementary energies that model both the spatial dependency of pixel classification in each follow-up exam and the information convoyed by all follow up scans (multisource fusion).Note that we have opted for the 8-connexity neighboring system in our MRF model of pixel classification.The proposed energy function is defined by: where p j,l (r(i, l)|q(i) = H j ) is the probability density function of the ith amplitude difference r(i) belonging to the source l (or HRT scan), where l = 1, 2, 3 and conditioned to H j ; α j,l is the set of the pdf hyperparameters; δ the delta Kroneker function; k = 1...N ; and γ l and β are positive parameters.Note that the first part of the energy function models the a priori we have on the image differences R(l) conditioned to change and no-change hypotheses which allows us to use the whole information available in every scan (information fusion).The second part models the spatial dependency by the use of the secondorder 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 classification results of the HRT images [24]s.The hyperparameters β, γ 1 , γ 2 , γ 3 handle the importance of the energy terms.On one hand, β allows us to tune the importance of the spatial dependency constraint.In particular, high values of β correspond to a high spatial constraint.On the other hand, γ l models the reliability of each HRT follow-up image and it is usually assumed to take on values in ]0, 1].This constraint can easily be satisfied by choosing the appropriate a priori distribution for γ l .Note that when γ l = 0, all the energy contribution related to the difference-images is removed and hence should be avoided.
Several parametric models can be used to model the distribution of r(i, l) conditioned to H j .In this work, we opted for the normal distribution.The pdf p j,l (r(i, l)|q(i) = H j ) is then given by: where µ j,l and σ j,l stand for the mean and the standard deviation respectively.
The whole set of the model hyperparameters is then given by Θ = {σ j,l , µ j,l , γ l , β}.In the absence of information prior knowledge, the following priors were used to generate model hyperparameters: The non-negativity of the the hyperparameters {γ l , β} is guaranteed through the use of the exponential densities: Note that the values of {η, ϕ, κ} are fixed empirically but do not really influence the results.
To estimate the model parameters and hyperparameters, we propose the use of a Monte Carlo Markov Chain (MCMC) procedure.Specifically, the principle of MCMC method is to generate samples drawn from the posterior densities and then to achieve parameter estimation.We use a Gibbs sampler based on a stationary ergodic Markov chain allowing us to draw samples whose distribution asymptotically follows the a posteriori densities.

C. Inference scheme
The Gibbs sampler facilitates efficient sampling from a n-dimensional joint distribution when knowledge about its conditional distributions is available.The aim of the approach is to generate a sample of values from the posterior distribution of the unknown parameters.The inference scheme consists of running the Gibbs sampler for many iterations and at each iteration an estimate for each unknown parameter is produced.At each iteration, the parameter estimate from the last iteration is used to produce new estimate by assuming that the current values for the other parameters are the true values.The main Gibbs sampler steps are described in algorithm 1.
j,l and µ [h] j,l from p(σ j,l ) and p(µ j,l ) iv) Sample γ For convergence, we used a burn-in period of h min = 500 iterations followed by another 1000 iterations for convergence (h max =1500).The change detection map Q is then estimated using the maximum a posteriori MAP estimator: Q = argmax Hj pHj , where pH0 = )) and pH1 = 1 − pH0 and p [h] (q i = H j |R, Θ)) is the estimated pdf at the iteration h.

III. EXPERIMENTS
This section aims at validating the proposed framework.Datasets used for validation are presented in the next subsection.In sub-section III-B, the intensity normalization algorithm reliability is presented.Change detection results on semi-simulated datasets are presented in sub-sections III-C.Finally, the glaucoma progression detection results using datasets are presented in sub-section III-D.

A. Datasets
The proposed framework was experimentally validated with clinical datasets.The three study groups in the clinical datasets have been described previously [6] .In brief, all eligible participants were recruited from the University of California, San Diego Diagnostic Innovations in Glaucoma Study (DIGS) with at least 4 good quality HRT-II exams, at least 5 good quality visual field exams and at least 2 good quality stereo-photographs of the optic disk were included in the study (267 eyes of 202 participants).Two hundred and forty six eyes from 167 glaucoma patients were included as progressing or non-progressing.Thirty six eyes from 33 participants progressed by stereo-photographs and / or showed likely visual field (Progressors) and the rest of the 210 eyes from 148 participants were considered nonprogressing (Non-progressors).An additional 21 eyes from 20 participants were normal eyes with no history of IOP>22 mmHg, normal appearing optic disk by stereo-photography and visual field exams within normal limits (median age of 62.7 years and median HRT-II follow-up of 0.5 years).The UCSD Institutional Review Boards approved the study methodologies and all methods adhered to the Declaration of Helsinki guidelines for research in human subjects and the Health Insurance Portability and Accountability Act (HIPAA).www.ijacsa.thesai.orgHRT exams from each study eye were of size 360 x 360 pixels.For each eye, several exams were performed over time.For each exam, the baseline image and the follow-up image are co-registered using built in instrument software and saved together.
Since no change detection map is available for the clinical dataset, we generated two semi-simulated datasets to assess the proposed change detection method.The first semi-simulated dataset (dataset1) was constructed from four normal HRT images.Changes were simulated by permuting 5%, 7%, 10% and 12% of image regions in the four images respectively (Figure 1).The second simulated dataset (dataset2) was constructed from another set of four normal HRT images.Changes were simulated by modifying the intensities of each 15 x 15 pixel-sized regions in the four images randomly with 0.5%, 0.75%, 1% and 1.25% probability of occurrence respectively.The intensities were modified by multiplying randomly the real intensities by 0.5, 0.75, 1.5 or 2 with a 0.25% probability (Figure 2).

B. Intensity normalization algorithm assessment
To assess the proposed illumination correction normalization algorithm and particularly the use of the retinex method for the reflectance estimation, five different methods were used to normalize the image intensities: the proposed method [18] (filter size=20 pixels), the homomorphic filtering method [16] (standard deviation=2 and filter size=20 pixels), the isotropic diffusion method [20] (smoothness constraint=7 pixels), the Discrete Cosine Transform (DCT) method [25] (number of components is 40) and the waveletbased method (with Daubechies wavelet and three level decomposition) [21].Fig. 3 presents the intensity normalization of a baseline and follow-up images using the five methods.
All the above methods were evaluated using the simulated datasets, dataset1 and dataset2.For evaluation, we use False Alarm PFA, Missed Detection PMD and Total Error PTE measurements computed in percentage and defined as: NM +NF × 100, where F A stands for the number of unchanged pixels that were incorrectly determined as changed, N F is the total number of unchanged pixels, M D the number of changed pixels that were incorrectly detected as unchanged, N M is the total number of changed pixels.Table .I presents false detections, missed detections and total errors for the simulated datasets.As one can see, the retinex reflectance algorithm performs better than the other normalization methods.

C. Change detection results
In order to emphasize the benefit of the proposed change detection algorithm and particularly the use of the Markov model to handle pixel spatial dependency, we compared the proposed method to the following two kernel-based methods and a Bayesian threshold-based method: • The Support Vector Data Description SVDD [26] with the Radial Basis Function RBF kernel, • The Support Vector Machine [27] with the RBF Gaussian kernel, • A Bayesian threshold-based method [28], Note that we used the retinex-based intensity normalization for all methods.First, we applied these methods on the semi-simulated datasets.The proposed method tends to perform better than the above methods.This means that the proposed Markov a priori we considered improves the change detection results by utilizing the information about spatial dependency of pixels in our detection scheme.

D. Glaucoma progression detection
We are now faced with the problem of framework validation on clinical datasets.We estimated the sensitivity and the specificity of detecting glaucoma progression as: sensitivity = T P T P +F N ; specif icity = T N T N +F P where T P stands for the number of true positive identifications, F N the number of false negative identifications and F P the number of false positive identifications.
In order to emphasize the benefit of the proposed glaucoma progression detection scheme, we have compared the proposed framework with three other published methods: the Topographic Change Analysis (TCA) method [5], the Proper Orthogonal Decomposition POD method [6] and the reflectance based method [7] (Table III).Note that an eye is considered as progressor when, in one exam, the number of changed pixels with an increase of intensity (loss of retinal height) is greater that 5% of the number of pixels within the optic disc [6].The proposed framework had approximately twice the specificity in the non-progressing eyes www.ijacsa.thesai.orgProgressor sensitivity (72%) than the POD method (43%).Moreover, it has higher sensitivity in progressor eyes (87%) while maintaining good specificity in normal eyes (91%).Further, in comparison to the reflectance based method, our proposed method provided similar specificity in normal and non-progressing eyes and higher sensitivity in progressing eyes.Increased sensitivity in our proposed scheme is likely because of the fact that we explicitly modeled the spatial dependency of clasification among pixels where as dependency is only implicitly accounted by the reflectance based method using familywise type I error rate.

IV. CONCLUSION
In this paper, a Bayesian framework for glaucoma progression detection has been proposed.The task of inferring the glaucomatous changes is tackled with a Monte Carlo Markov Chain algorithm that is used for the first time to our knowledge in the glaucoma diagnosis framework.Modeling and accounting of both spatial dependency of pixels increased the robustness of the proposed change detection scheme compared to the kernel-based method and the threshold method.The validation of the proposed approach using clinical datasets has shown its ability to provide high specificity in non-progressor stable glaucoma eyes and high sensitivity in progressor eyes while maintaining a good specificity in normal eyes.

Fig. 1 .Fig. 2 .
Fig. 1.Semi-simulated image by the permutation of image regions: (a) the original image (b) the simulated image (PSNR=28 dB) and (c) the ground truth.

Fig. 3 .
Fig. 3. Intensity normalization results of a baseline and a follow-up images respectively using : (a) (b) the proposed method, (c) (d) the homomorphic filtering, (e) (f) the isotropic diffusion method, (g) (h) the DCT method and (i) (j) the wavelet method.228| P a g e www.ijacsa.thesai.org