Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2749882

Formats

Article sections

- Abstract
- I. Introduction
- II. Background
- III. ONH Topographic Library from the LSU Experimental Glaucoma Study
- IV. Rotational Alignment of TopSS Topographies
- V. Building an Optimal Baseline-subspace of the ONH of an Eye Using POD
- VI. Glaucomatous Change Detection Framework
- VII. Performance Analysis
- VIII. Results
- IX. Discussion
- X. Conclusion
- References

Authors

Related links

IEEE Trans Inf Technol Biomed. Author manuscript; available in PMC 2010 August 1.

Published in final edited form as:

Published online 2009 April 14. doi: 10.1109/TITB.2009.2020158

PMCID: PMC2749882

NIHMSID: NIHMS110296

Madhusudhanan Balasubramanian, Member, IEEE,^{1} Stanislav Žabić,^{2} Christopher Bowd,^{1} Hilary W. Thompson,^{3} Peter Wolenski,^{2} S. Sitharama Iyengar, Fellow, IEEE,^{4} Bijaya B. Karki,^{4} and Linda M. Zangwill^{1}

Madhusudhanan Balasubramanian: ude.dscu.amocualg@uhdam; Stanislav Žabić: moc.liamg@1cibazs; Christopher Bowd: ude.dscu.amocualg@dwobc; Hilary W. Thompson: ude.cshusl@2pmohTH; Peter Wolenski: ude.usl.htam@iksnelow; S. Sitharama Iyengar: ude.usl.csc@ragneyi; Bijaya B. Karki: ude.usl.csc@ikrak; Linda M. Zangwill: ude.dscu.amocualg@lliwgnaz

The publisher's final edited version of this article is available at IEEE Trans Inf Technol Biomed

See other articles in PMC that cite the published article.

Glaucoma is the second leading cause of blindness worldwide. Often the optic nerve head (ONH) glaucomatous damage and ONH changes occur prior to visual field loss and are observable *in vivo*. Thus, digital image analysis is a promising choice for detecting the onset and/or progression of glaucoma. In this work, we present a new framework for detecting glaucomatous changes in the ONH of an eye using the method of proper orthogonal decomposition (POD). A baseline topograph subspace was constructed for each eye to describe the structure of the ONH of the eye at a reference/baseline condition using POD. Any glaucomatous changes in the ONH of the eye present during a follow-up exam were estimated by comparing the follow-up ONH topography with its baseline topograph subspace representation. Image correspondence measures of *L*_{1} and *L*_{2} norms, correlation, and image Euclidean distance (IMED) were used to quantify the ONH changes. An ONH topographic library built from the Louisiana State University Experimental Glaucoma study was used to evaluate the performance of the proposed method. The area under the receiver operating characteristic curves (AUC) were used to compare the diagnostic performance of the POD induced parameters with the parameters of *Topographic Change Analysis* (TCA) method. The IMED and *L*_{2} norm parameters in the POD framework provided the highest AUC of 0.94 at 10° field of imaging and 0.91 at 15° field of imaging compared to the TCA parameters with an AUC of 0.86 and 0.88 respectively. The proposed POD framework captures the instrument measurement variability and inherent structure variability and shows promise for improving our ability to detect glaucomatous change over time in glaucoma management.

Glaucoma is a progressive optic neuropathy that, when left untreated, may result in progressive vision impairment and eventual blindness [1]. It is the second leading cause of blindness in the world next only to cataract [2]. Current estimates based on population based surveys predict that there will be 60 million people worldwide affected by glaucoma in 2010, increasing to 80 million by 2020 [3]. Retrospective estimates of glaucoma related expenses such as ophthalmologist visits, glaucoma related surgical procedures, medications, and other indirect costs and services indicate an annual cost of about $1.2 billion to $2.5 billion in the US [4]–[6].

The retina is a photoreceptive layer that is composed, in part, of about one million optic nerve fibers originating from the neuron cell bodies within the ganglion nerve fiber layer. The nerve fibers group together and form the nerve fiber bundle or optic nerves. These optic nerve fibers run parallel to the retinal surface and exit the eye at the optic disk (blind spot) located in the posterior end of the eye and carry visual electrical impulses to the visual cortex of the brain for image formation. The optic disk region is usually referred to as the optic nerve head (ONH) and exhibits a natural cup shape due to the arrangement of the optic nerves leaving the eye. The natural arrangement of the optic nerves as they exit the eye results in a sloping region called the neuroretinal rim and a deeper region called as an optic cup. Fig. 1 shows the ONH topography of a primate subject from the Louisiana State University Experimental Glaucoma study.

ONH topographies of subject 1D indicating differences in the resolution at 10° and 15° field of imaging. The topographic surface represents the vitreo-retinal interface of the eye. An ellipse is manually fit to mark the optic disk margin **...**

A detailed background on the pathophysiology of glaucoma is described elsewhere [7]. Because glaucoma is a chronic disease, it typically results in a gradual loss of nerve fibers associated with retinal ganglion cells. The loss of nerve fibers causes characteristic changes in the appearance of the retinal nerve fiber layer and eventual changes in the configuration of the optic disk. Therefore, digital image analysis of the optic disk region is a rapid and promising approach for detecting glaucomatous changes in the ONH region of an eye. Confocal Scanning Laser Ophthalmoscope (CSLO), a class of confocal microscope that utilizes a rotating mirror arrangement to scan an imaging area at various z-axis depths using laser light beams, can be used to capture the 3-D architecture of the ONH region. The Heidelberg Retina Tomograph (HRT; Heidelberg Engineering, GmbH, Heildelberg, Germany) and Topographic Scanning System (TopSS; formerly of Laser Diagnostic Technologies, San Diego, CA) are two of the CSLO instruments used for imaging the ONH region.

In the early stages of glaucoma, the characteristic changes associated with structural glaucomatous damages may not be drastic and therefore, may not be obvious. Because glaucoma is characterized as a progressive optic neuropathy, the diagnostic accuracy can be improved by detecting changes in the ONH structure of an eye from a reference or baseline condition. For detecting progressive glaucomatous changes, the ONH topographies constructed from the follow-up CSLO exams are compared with the baseline topographies of the eye. Two methods currently available for detecting topographic pixel-level glaucomatous changes are: 1) a nonparametric permutation test based method known as the Statistic Image Mapping (SIM) of the retina [8], and 2) an ANOVA model based change detection method [9] available in the HRT software as Topographic Change Analysis (TCA).

In this work, we present a novel computational framework for detecting glaucomatous changes in the ONH region of an eye at the original topographic resolution. The method was inspired from the theory of proper orthogonal decomposition (POD) and its application in the analysis of dynamics of turbulence in fluid mechanics [10]. POD is theoretically similar to Karhunen-Loève expansion, Hotélling transform, and principal component analysis (PCA). Previous related use of these techniques include facial recognition application [11]–[13], video scene change detection [14], estimating brain volume changes in patients with probable Alzheimer’s disease [15], restoration of white-light confocal microscope optical section images [16] and for detecting changes in the satellite images [17].

A brief background of the SIM and TCA methods is presented in Section II. The proposed POD framework was tested using a library of ONH topographies of 12 primates built from the Louisiana State University Experimental Glaucoma (LEG) study [18]. A brief description of the ONH topographic library from the LEG study is presented in Section III. In Section IV, we describe a semi-automatic optic disk rotational alignment algorithm using an FFT based template matching algorithm that we used to correct for any rotational alignment errors in the ONH topographies. Details of the POD glaucomatous change detection framework are presented in Sections V and VI. The performance of the POD framework is compared with the TCA and the results are presented in Section VIII. For demonstration, we use the ONH topographies of subject 1D from the LEG study in all the figures (except Fig. 3).

The SIM of the retina method is based on suprathreshold cluster tests [19] wherein the locations with statistically significant changes are chosen using a primary threshold (for example, the locations with *p* < 0.05). Because glaucomatous structural changes affect a contiguous region in the ONH, the locations with significant changes are further grouped into clusters and the significance of the topograph-level progression are defined based on the significance of the size of the largest cluster of significantly changed locations using a non-parametric permutation test. The SIM of the retina method utilizes a slope-based test-statistic to estimate the primary pixel-level changes in the topographic height measurements using a baseline and several follow-up visits. Significance of the pixel-level test-statistics and the significance of the largest cluster of significantly changed locations are estimated using nonparametric permutation tests. The current SIM of the retina method requires a minimum of 7 follow-up exams to build a reliable permutation distribution for the slope test-statistic and for the cluster-size significance test.

The HRT TCA method utilizes a mixed-effect three-way ANOVA model to detect *superpixel* level changes in the ONH topographies from a baseline exam to a follow-up exam. A superpixel is a group of topographic height measurements from a neighborhood of 4 ×4 individual pixels pooled in the ANOVA model for a specific detection of localized changes. The ANOVA model for detecting superpixel changes, with each superpixel containing 16 topographic height measurements from locations *l* = 1, 2, …, 16, from a baseline (at time *t* = 1) to a follow-up visit (at time *t* = 2) using a set of N topographies each acquired at the baseline and follow-up visits, is given by

$${h}_{\mathit{tli}}=\mu \dots +{T}_{t}+{L}_{l}+I{(T)}_{i(t)}+T{L}_{tl}+{\epsilon}_{\mathit{tli}}$$

(1)

where, *h _{tli}* is the retinal height at location

The significance of the mean retinal height change at each superpixel location is estimated using Satterthwaite’s approximate F-test. For each follow-up exam, a change significance map is constructed by identifying the superpixel locations with significant decrease in retinal height from the baseline exam (i.e., locations with negative height change in the mean difference topography and with change probability *p* < 0.05). Fig. 2 demonstrates the application of the TCA using the TopSS topographies of a primate eye with experimentally induced glaucoma from the LEG study.

Detecting significant retinal height changes from the baseline session N1 to the follow-up session 01 of subject 1D using topographic change analysis (TCA)

For a specific detection of glaucomatous changes in a follow-up exam, superpixel locations with significant decrease in retinal height repeatable in 2- or 3-successive follow-up exams are grouped into clusters [20]. In TCA, commonly used summary parameters estimated based on changes repeatable in 2- or 3-successive follow-up exams are: 1) size of the largest cluster of superpixels with significant height decrease within the optic disk margin (CSIZE), 2) proportion of the CSIZE measure to the optic disk size (CSIZE%), and 3) total number of superpixels with significant height decrease within the optic disk margin [20]–[22].

Glaucomatous changes in the ONH region, that may result in changes in the optic disk configuration and/or in the nerve fiber layer, originate either directly from the defects in the underlying neuronal cell layers or indirectly influence the function of the neural cell layers and cause eventual changes in the visual function of the eye. However, the order of the temporal sequence of observing the ONH structural changes and the visual function defects in an eye vary [23]. Although the visual function defects are associated with structural defects such as dysfunction and death of retinal ganglion cells and loss of optic nerve fibers, the associated structural defects may not always be visible in the retinal surface at the same time; i.e. the topographic ONH surface changes are not temporally correlated with the visual function changes. For example, when the initial damages occur in the deeper underlying layers of the retina, changes in the visual function of the eye may be observed earlier than the appearance of the structural changes in the ONH topographies. Due to this difficulty, currently, there is no single gold standard available to define a glaucomatous progression in an eye. Current studies utilize an operator evaluation of optic disk photographs and/or changes in the visual function to define glaucomatous progression in an eye. Therefore, to evaluate a new glaucomatous progression algorithm, it is an ideal choice to experimentally induce glaucoma in an eye, wherein the state of an ONH is experimentally controlled during each follow-up exam. In experimental glaucoma, the IOP of an eye is elevated above the normal level to induce optic nerve damages.

In this study, we utilize a library of ONH topographies of 12 primates (24 eyes), built from the LEG study [18], to evaluate the performance of the POD framework. Details of all aspects of the LEG study have been described previously [18]. In brief, one eye of each primate was treated with laser to induce glaucoma (*glaucoma induced study eye*) and the other eye is untreated (*contralateral normal eye*). The ONH of both the glaucoma induced study eye and the contralateral normal eyes of all the primates were imaged every two weeks. During each imaging sessions, 6 ONH topographies of each eye were obtained to characterize the topographic measurement variability of an eye within each imaging session. The IOP in the glaucoma induced study eye of each of the 12 primates were elevated by treating the trabecular meshwork with argon laser. Prior to inducing IOP elevation in the glaucoma induced study eyes, the ONH of both the eyes were imaged during three separate baseline imaging sessions (*N*_{1}, *N*_{2}, *N*_{3}) spaced at least two weeks apart. After elevating the IOP level in the glaucoma induced study eyes, both the glaucoma induced study eyes and contralateral normal eyes were imaged every two weeks as mentioned earlier. A summary of the pre- and post-laser imaging sessions of all the primates at 10° and 15° field of imaging is presented in Table I. Out of the 6 topographic scans acquired during each imaging session, a minimum of 4 topographies were required to be of good quality and the poor quality scans were discarded. Imaging sessions with less than 4 good quality scans were repeated. A total of 2,098 good quality ONH topographies at 10° field of imaging and 2,107 good quality ONH topographies at 15° field of imaging are available in the LEG study library.

TopSS—a confocal scanning laser ophthalmoscope (CSLO), was used to image the ONH of the primates. Principle of a CSLO and the details of TopSS have been described previously [18]. In brief, TopSS has a lateral digital resolution ranging from 11*μ*m at 10° field of imaging to 35*μ*m at 30° field of imaging and an axial resolution of ~300*μ*m (full width at half maximum). Reproducibility of topographic measurements are 20*μ*m laterally and better than 30*μ*m axially. TopSS constructs a topography that represents the surface of the ONH from a set of 32 optical serial section images of the ONH region. The optical section images, each of size 256 ×256 pixels, correspond to the ONH structure at various depths and collectively represent the 3D structure of the ONH. A least squares fit algorithm is used to register the optical section images of a scan using one of the optical sections from the scan as a reference [18]. An ONH topography is constructed by identifying the maximal reflectance value at each pixel position among all of the 32 registered optical section images (see Fig. 3 for an example of constructing a ONH topography from optical section images). Thus, a topographic surface is expected to represent the vitreo-retinal interface in an eye. Topographies acquired during follow-up imaging sessions are registered to one of the baseline topographies using a least squares fit algorithm [18]. The imaging angle of the instrument can be set to acquire images at either 10 ° or a 15 ° field of imaging while maintaining the same digital imaging area of 256 ×256 pixels. ONH topographies at 10 ° field of imaging are at a higher resolution compared to the topographies at 15 ° field of imaging and therefore, incorporate more ONH details as shown in Fig. 1.

The TopSS software allows studying changes in the optic disk configuration by manually drawing an ellipse on the outer margin of the neuroretinal rim or the inner margin of the scleral ring in the baseline exam of an eye. In the LEG study, a TopSS operator manually chose one of the ONH topographies from the baseline imaging session of an eye as a reference and selected ellipse parameters *a* and *b* and coordinates (*x*_{0}, *y*_{0}) appropriate for the ONH being imaged to manually mark the optic disk margin. Because the shape of an optic disk is vertically oval [7], TopSS aligns the major axis of the ellipse along the y-axis (i.e. at an orientation angle of 90°). The ellipse parameters were automatically transferred to each of the follow-up ONH topographies upon registration and is useful in determining the glaucomatous changes over time in the neuroretinal rim and optic cup. Further, the operator selected a set of new ellipse parameters and coordinates when there was a significant change in the shape of the optic disk during a follow-up or when the ellipse transferred to a follow-up topography by the TopSS software did not adequately fit the optic disk margin. Ellipse parameters from the most recent imaging session of an eye, approved by the TopSS operator, were applied to the subsequent follow-up exams of the eye.

Our evaluation of the time series of TopSS ONH topographies from the LEG study library revealed rotational misalignment errors as large as ~18° among baseline and follow-up topographies of an eye; for example see Figs. 4a, 4c and 4f for the topographies of subject 1D from imaging sessions *N*_{1}, *N*_{2} and 06 respectively. The operator evaluation of the ellipses transferred to the follow-up topographies guaranteed a more accurate ellipse fit for the optic disk and more importantly allowed to correct any translational errors (i.e. vertical and horizontal shifts) among ONH topographies after automated alignment of topographies using the TopSS software. However, the topographies from the library revealed that, this procedure did not adequately account for any rotational errors not corrected by the TopSS registration algorithm. Because both the POD framework and the TCA method detect localized pixel-level changes, their detection accuracies are dependent on the alignment of the baseline and follow-up ONH topographies of an eye. Therefore, we corrected any rotational misalignment among the ONH topographies using an FFT based template matching algorithm and the details are presented in the next section.

Rotational misalignment between a baseline and a follow-up topography with a common center, can be estimated as a translation along the *θ*-axis in the polar coordinates [24]. Topographies, in the cartesian coordinate system, can be transformed to polar coordinates by interpolating radial pixel intensities along the *r*-axis at various angles *θ* from the topographic height matrix. For example, pixel intensity at a radius *r* and an angle *θ* can be determined by interpolating the pixel intensity at (*x* = *r* cos *θ*, *y* = *r* sin *θ*) from the topography. We used bilinear interpolation for transforming ONH topographies to polar coordinates. Fig. 4a shows an ONH topography of subject 1D from baseline session *N*_{1} (in cartesian coordinate) and Fig. 4b shows its corresponding polar coordinate representation for a circular region shown in Fig. 4a. Manual evaluation of the ellipse fit for the optic disk margin by a trained TopSS operator, as mentioned in section III, corrects for any horizontal and vertical shifts of the optic disk region in the follow-up imaging sessions with reference to a baseline topography.

The optic disk region (ODR) enclosing the profiles of peripapillary arteries and veins in ONH topographies were used to estimate the rotational misalignment errors. The rotational misalignment angle in a follow-up ODR can be determined using a *template matching* algorithm [25]. Let *I _{bl}* and

$$\frac{{\sum}_{j=1}^{Q}{\sum}_{k=1}^{R}{I}_{t}(j,k){\widehat{I}}_{\mathit{fup}}({i}_{r}+j,{i}_{\theta}+k)}{\sqrt{{\sum}_{j=1}^{Q}{\sum}_{k=1}^{R}{I}_{t}^{2}(j,k)}\sqrt{{\sum}_{j=1}^{Q}{\sum}_{k=1}^{R}{\widehat{I}}_{\mathit{fup}}^{2}({i}_{r}+j,{i}_{\theta}+k)}}$$

(2)

where, *i _{r}* and

A faster, yet, a reliable template matching method using a normalized correlation measure as in (2) can be implemented using FFT [25], [26]. With representing the Fourier transform, the numerator of (2) can be efficiently implemented as a convolution of *Î _{f up}* with

$$\frac{{\mathcal{F}}^{-1}\left(\mathcal{F}({I}_{t})\mathcal{F}{({\widehat{I}}_{\mathit{fup}})}^{\ast}\right)}{\sqrt{{\sum}_{j=1}^{M}{\sum}_{k=1}^{N}{I}_{t}^{2}(i,j)}\sqrt{{\mathcal{F}}^{-1}\left(\mathcal{F}({\widehat{I}}_{\mathit{fup}}^{2})\mathcal{F}{(\mathit{Mask})}^{\ast}\right)}}$$

(3)

Our implementation of the FFT based normalized template matching algorithm in MATLAB ver. 7.6 provided a significant improvement in computational speed compared to the direct normalized correlation based template matching technique. For e.g., searching a template of size 50×50 pixels in an image of size 256×256 pixels took 0.056 seconds using the FFT based method compared to 4.8 seconds using the direct method.

For estimating true translations along the *θ*-axis (i.e. rotational error), the template chosen for matching should enclose spatial relationships between features that are unique and stable during follow-up imaging sessions. The spatial arrangement of arteries and veins branching out from the optic disk region is relatively stable during the disease process and therefore, can be utilized to identify rotational misalignment errors. In contrast, due to glaucomatous progressive changes during follow-ups, the shape of the optic disk may change significantly from a baseline condition. Therefore, to improve the rotational error estimates, we excluded the optic disk region and utilized arteries and veins in the peripapillary retina while choosing a template *I _{t}* from

A proper orthogonal decomposition of an ensemble of vectors provides an optimal set of orthogonal basis vectors that span the entire ensemble space [27] and guarantees the best *k*-term orthogonal expansion, in a mean squared error sense, among all the orthogonal transforms [28]. POD decorrelates a given ensemble of signal by discovering an orthogonal basis set that is specific and optimal for the signal under consideration as opposed to analyzing the signal using an off-the-shelf wavelet basis or Fourier transform.

Let,
${\{{I}_{bl}(m)\}}_{m=1}^{N}$be a set of expectation-centered grayscale ODRs from the ONH topographies acquired at a baseline or reference condition of an eye *E _{s}* of a given subject

Using POD, each of the ODRs *I _{bl}*(

$$\begin{array}{l}{I}_{bl}(m)=\sum {\alpha}_{(m,n)}{\phi}_{n}\\ \text{with},\phantom{\rule{0.38889em}{0ex}}E({\alpha}_{(m,i)},{\alpha}_{(m,j)})=\sqrt{{\lambda}_{i}}{\delta}_{ij}\phantom{\rule{0.16667em}{0ex}}\text{and}\\ \langle {\phi}_{i},{\phi}_{j}\rangle ={\delta}_{ij}\end{array}$$

(4)

Here, *E* is the expectation, {*λ _{i}*} are the eigenvalues and {

$$R\mathrm{\Phi}=\mathrm{\Lambda}\mathrm{\Phi}$$

(5)

Thus, the eigenvectors Φ computed using the covariance *R* of the ensemble matrix *X* form an optimal orthogonal basis for the baseline ODR ensemble. The existence of Λ = {*λ _{i}*} and Φ = {

The eigenvectors Φ = {* _{i}*} in (5) that form the orthogonal basis for the ensemble can be derived using a singular value decomposition (SVD) of the ensemble matrix

$$X=U\mathrm{\Lambda}V$$

(6)

$$\text{Covariance}\phantom{\rule{0.16667em}{0ex}}R=X{X}^{T}=U\mathrm{\Lambda}V{V}^{T}\mathrm{\Lambda}{U}^{T}=U{\mathrm{\Lambda}}^{2}{U}^{T}$$

(7)

The left singular vectors *U* span the column of the ensemble matrix *X* and thus, form the basis of the ODR ensemble. Solving the eigenvalue problem in (7), for the ODR ensemble matrix *X* formed using the baseline ensemble {*I _{bl}*(

$$\begin{array}{l}{X}^{T}X={V}^{T}\mathrm{\Lambda}{U}^{T}U\mathrm{\Lambda}V\\ ={V}^{T}{\mathrm{\Lambda}}^{2}V\end{array}$$

(8)

The right singular vectors *V* and the eigenvalues Λ can be computed from (8). This requires solving a *N* ×*N* system. Further, the left singular vectors *U* can be computed from (6) using *X*, *V* and Λ. The reduced SVD approach for determining an optimal ensemble basis provides a significant reduction in computation when the ensemble size *N* *ST*, with an ODR size of *S* ×*T*.

Prior to applying the change detection framework described in this section, rotational misalignments that we observed in the ONH topographic library from the LEG study were corrected using the rotational alignment procedure described in section IV. In general, after delineating the ODR of an eye by manually drawing a contour line (contour line is an ellipse in the case of TopSS topographies and a spline curve in the case of HRT topographies) and after aligning baseline and follow-up topographies of an eye, the following steps can be applied to ONH topographies obtained at either 10° or 15° field of imaging to identify changes in the ODR of the eye from baseline.

We selected ODR from each of the topographies by constructing a minimum bounding rectangle of size (2*a* + 1) × (2*b*+1) with center (*x*_{0}, *y*_{0}) guided by the ellipse coordinates (*x*_{0}, *y*_{0}) and parameters *a* and *b* that mark the optic disk margin in a topography. To select the same area of ODR from all the exams of an eye for analysis, we chose ellipse parameters *a* and *b* from one of the baseline exams of the eye and applied to each of the follow-up exams. Alignment among baseline and follow-up exams were maintained by using their respective ellipse center coordinates (*x*_{0}, *y*_{0}). Because any rotational alignment errors were corrected before this step, using the same baseline ellipse parameters *a* and *b* for all the follow-up imaging sessions ensured that a similar ODR was selected for analysis from all the imaging sessions of an eye.

As mentioned earlier, let *I _{bl}* and

An optimal orthogonal subspace at baseline was constructed for each eye as described in section V using the ODRs of all the baseline topographies of the respective eye (ODR selected as in step *A*). The optimal orthogonal subspace, called *baseline subspace*, describes the baseline condition of an eye and is used to quantify changes in the ONH of the same eye observed during a follow-up at a later time. The baseline topographies chosen for subspace construction should characterize the variability in the appearance of the ONH at the baseline condition and the measurement variability due to the imaging instrument. In the LEG study, the baseline imaging sessions *N*_{1}, *N*_{2} and *N*_{3} were scheduled at least two weeks apart before inducing an experimental glaucoma in the glaucoma induced study eyes in order to capture the variability in the ONH structure at the baseline condition. Therefore, we used topographies from all the 3 baseline imaging sessions of an eye to construct the baseline subspace of the respective eye.

Let
${\{{\phi}_{i}\}}_{i=1}^{N}$represent the baseline subspace of an eye, where *N* is the size of the ensemble of baseline topographies used for subspace construction. Fig. 5a shows an ensemble of baseline ONH topographies from the pre-laser treatment imaging sessions *N*1, *N*_{2}, and *N*_{3} of subject 1D and Fig. 5b shows a pictorial representation of a set of optimal POD basis vectors that form the baseline subspace of the eye.

To determine changes in the ODR of a follow-up topography *I _{f up}*, we identified a topography

$${\widehat{I}}_{\mathit{fup}}=\sum _{i=1}^{N}\langle {I}_{\mathit{fup}},{\phi}_{i}\rangle {\phi}_{i}$$

(9)

Fig. 6a shows the mean topographies from the pre- and post-laser treatment imaging sessions of the contralateral normal eye of subject 1D and Fig. 6b shows their respective baseline subspace representations. Fig. 7a and 7b show the mean topographies and their baseline subspace representations for the glaucoma induced study eyes of subject 1D. The corresponding TCA change significance maps for the contralateral normal eye and the glaucoma induced study eye of subject 1D are shown in Figs. 6c and and7c7c respectively. For the contralateral normal eye, it can be seen that the POD baseline subspace representation of the ONH topographies from the follow-up imaging sessions (Fig. 6b) closely represent the respective original ONH topographies (Fig. 6a). Because the laser-treated eyes changed significantly due to the experimental glaucoma from the pre-laser treatment baseline condition, the POD baseline subspace cannot accurately describe their post-laser treatment ONH topographies as shown in Fig. 7b.

POD analysis and TCA of the contralateral normal eye of Subject 1D. The observed topographies in Fig. a and their respective baseline subspace representations in Fig. b appear more similar indicating less changes from baseline. For e.g., IMED 6,338 of **...**

It can be observed that when there are minimal or no changes in a follow-up exam from a baseline exam, the baseline subspace of the eye can accurately describe the corresponding follow-up topography. Therefore, glaucomatous changes in the follow-up ODR *I _{f up}* were quantified by determining the correspondence in topographic measurements from locations with decrease in retinal height between follow-up

*L*_{1}norm,*Norm*_{L1}(*I*,_{f up}*Î*), computed as_{f up}$$\sum _{i}\sum _{j}\mid {I}_{\mathit{fup}}(i,j)-{\widehat{I}}_{\mathit{fup}}(i,j)\mid $$(10)*L*_{2}norm or Euclidean distance,*Norm*_{L2}(*I*,_{f up}*Î*), computed as_{f up}$$\sqrt{\sum _{i}\sum _{j}{({I}_{\mathit{fup}}(i,j)-{\widehat{I}}_{\mathit{fup}}(i,j))}^{2}}$$(11)- Normalized correlation coefficient computed as in (2)

The *L*_{1} and *L*_{2} norms give a measure of mismatch between *I _{f up}* and

Performance degradation in using Euclidean distance for image similarity measurement comes primarily from the pixel-to-pixel correspondence used and is due to the orthogonal coordinate system employed for measuring an image distance. The IMED overcomes this drawback by assigning a varying weight to the adjacent pixels using a non-orthogonal basis. To account for small pixel displacements, the IMED uses a nonorthogonal basis that assigns a varying weight to the adjacent pixels while comparing a pixel value between images. Because insignificant pixel movements are expected to be near their respective original pixel locations, a Gaussian kernel based pixel weighing scheme would be a natural choice for the image similarity measurement. The problem of determining an optimal non-orthogonal basis for computing the IMED is avoided by using a metric coefficient matrix *G* induced from the non-orthogonal basis. Using the symmetric positive definite matrix *G*, the IMED can be computed as a *G*-innerproduct as follows.

$$\begin{array}{l}\text{IMED}({I}_{\mathit{fup}},{\widehat{I}}_{\mathit{fup}})={\langle {I}_{\mathit{fup}},{\widehat{I}}_{\mathit{fup}}\rangle}_{G}\\ ={I}_{\mathit{fup}}G{\widehat{I}}_{\mathit{fup}}\end{array}$$

If the standard deviation *σ* of the Gaussian kernel is far less than the dimension of *I _{f up}*, then the construction of the metric coefficient matrix

$$\text{IMED}({I}_{\mathit{fup}},{\widehat{I}}_{\mathit{fup}})=\langle {I}_{\mathit{fup}},{\stackrel{\sim}{I}}_{\mathit{fup}}\rangle $$

(12)

where *Ĩ _{f up}* can be computed by filtering

For evaluating the performance of the POD framework and TCA, we utilized estimates of area under their respective receiver operating characteristic curves (AUC). The AUC gives a measure of average sensitivity of a diagnostic test for all specificity values and thus, provides an overall diagnostic performance of the test. Refer to [32] for a recent literature survey on receiver operating characteristic (ROC) curves and associated diagnostic measures.

For constructing the ROC curves of each of the progression summary parameters in the POD framework and TCA, the respective progression summary parameters from each of the post-laser treatment imaging sessions were grouped into a stable group and a progressing group. The post-laser treatment follow-up imaging sessions (i.e. imaging sessions 01 and onwards) from the glaucoma induced study eyes of all the primates were considered to be progressed and the follow-ups from the contralateral normal eyes without any laser treatment were considered to be stable. Using the DeLong’s non-parametric method of comparing correlated ROC curves, we estimated differences in the AUC estimates (and 95% confidence intervals) between the best performing summary parameter in the POD framework and one of the best performing summary parameters in the TCA method. [33], [34]

Fig. 9 shows the class-conditional nonparametric probability densities of the POD IMED parameter and the TCA CSIZE parameter for the stable and progressing group of ONH imaging sessions at 10° and 15° field of imaging. The probability densities for the progression summary parameters were built using the adaptive mixtures method [35].

Stable and progression class-conditional probability density functions of the POD IMED parameter and the TCA CSIZE parameter.

The ROC curves of the POD IMED parameter and the TCA CSIZE parameter at 10° and 15° field of imaging are shown in Fig. 10. The AUC of all the POD parameters and TCA parameters are listed in Table II. The IMED and *L*_{2} norm parameters in the POD framework resulted in the highest AUC of 0.94 at 10° field of imaging and 0.91 at 15° field of imaging among all of the progression summary measures. At 10° field of imaging, the difference in the AUC of 0.08 (95% CI = (0.04, 0.11); *p*-value < 0.0001) between the IMED/*L*_{2} norm parameters and the best performing TCA CSIZE/CSIZE% parameters was statistically significant. The respective AUC difference of 0.03 (95% CI = (0.00, 0.06); *p*-value = 0.063) at 15° field of imaging was of borderline significance.

Comparative ROC curves and AUC estimates of the POD IMED parameter and the TCA CSIZE parameter; ROC curves of other parameters were not plotted for clarity

Glaucoma is a disease of progression characterized by progressive changes in the ONH structure and/or visual function of the eye. Glaucomatous structural changes are often characteristic in the neuroretinal rim, optic cup and the nerve fiber layer, while the exact role of relevant structures such as *lamina cribrosa* are still under investigation. The POD framework presented in this study focusses on detecting changes exhibited in the optic disk region in the neuroretinal rim and optic cup and provides the highest diagnostic AUC measures, statistically significant at 10° field of imaging, for the experimental glaucoma primate subjects from the LEG study. Also the class-conditional probability densities of the POD IMED parameter has the least error probability of 0.14 (estimated as the area of overlap between the probability density functions of the stable and progression groups) at 10° field of imaging. The least error probabilities associated with the POD parameters show promise for further performance improvements using the POD framework at 10° field of imaging. Error probabilities at 15° field of imaging are generally higher with a decrease in the AUC of 0.91 from 0.94 at 10° field of imaging. This reduction in AUC and the associated increase in error probability is likely due to the lower imaging resolution available at the 15° field of imaging and therefore, it is preferable to acquire and use high resolution ONH topographies for analysis using the POD framework.

The pixel-wise change analysis technique of TCA analyze ONH topographies at 1/4th of the original topographic resolution. Also TCA requires two or more additional follow-up exams to confirm and establish confidence on a detected change. In contrast, the POD framework analyzes the ONH structure at the original topographic resolution and detects glaucomatous changes from a baseline condition using only a single follow-up exam.

Both the POD framework and the TCA method are sensitive to any topographic misalignments because the analysis is carried out at the pixel level. The IMED parameter in the POD framework can account for small topographic misalignments. In contrast to the TopSS instrument, the latest HRT-3 instruments acquire high resolution topographic scans for analysis (for e.g., compare TopSS topographies in Fig. 1 with HRT topographies in Fig. 3c) and utilize an improved feature based alignment algorithm. With a more accurately aligned high resolution ONH topographies from HRT, we expect a similar or an improved diagnostic performance of the POD framework. Performance of the proposed POD framework in a larger group of human participants from a population based glaucoma study will be investigated in a separate work.

Due to the anatomical arrangement of the ganglion cell axons exiting the eye, neuroretinal rim is typically thicker in the inferior, superior, nasal and temporal region within the optic disk, in that order, thus giving a characteristic shape to the optic disk and optic cup. Therefore, detecting changes within the optic disk margin is expected to provide a specific detection of glaucomatous changes in an eye. Other techniques such as TCA and SIM of the retina methods define glaucomatous changes based on observed pixel changes within the optic disk margin. Therefore in this work, we derived glaucomatous change summary parameters in the POD framework based on topographic measurements within the optic disk margin. However, in the POD framework, similar change summary measures can also be estimated in the peripapillary retina for detecting nerve fiber layer defects.

In the current analysis, the POD framework does not have a graphical representation of pixel-wise change locations. However, one of the advantages of the POD framework is that other statistical and computational pixel-wise change detection algorithms can be integrated with the POD framework. For example, after constructing the baseline subspace representations of follow-up topographies in section VI-C, pixel-wise changes between follow-up topographies and their respective baseline subspace representations can be estimated using a statistical procedure as in the TCA method. Inducing other pixel-wise change detection algorithms within the POD framework will be studied separately in a future work.

In TopSS, the optic disk region of an eye is delineated using a manually drawn ellipse with its major axis along the y-axis. When the optic disk is tilted and not aligned along y-axis, the ellipse fit may not accurately delineate the optic disk region of an eye. However, the effect of this limitation on the current analysis is minimal because a similar region was chosen from each of the baseline and follow-up exams of an eye for detecting change over time. HRT instruments overcome this limitation by using flexible spline curves to accurately delineate the optic disk region. Using the spline curve coordinates, optic disk region from HRT topographies can be selected for analysis by constructing a minimum bounding rectangular region as described in section VI-A.

We have presented a novel subspace approach for detecting glaucomatous progression in the ONH region of an eye. An advantage of the POD framework is its ability to characterize the topographic measurement variability due to the imaging instrument, and imaging conditions and more importantly characterize the inherent variability of an ONH structure at a reference or baseline condition. The POD framework constructs an optimal baseline subspace from a set of baseline topographies to characterize the baseline variability of an eye. Glaucomatous changes in a follow-up exam, typically beyond the structural and measurement variability observed at baseline, can be detected by comparing the follow-up topography with the baseline subspace. The IMED and *L*_{2} norm parameters in the POD framework provide the highest diagnostic accuracies at both 10° and 15° field of imaging. The POD framework provides the best diagnostic accuracy when using high resolution topographies (i.e. TopSS scans at 10° field of imaging) and therefore it is preferable to use high resolution scans for analysis using the POD framework. The proposed POD framework can also be used with other imaging modalities such as optical coherence tomography and scanning laser polarimetry for detecting glaucomatous structural change over time and can also be extended for detecting visual function changes.

This research was supported in part by the Grant Number P20RR016456 from the National Center For Research Resources, the Grant Number EY11008 (LMZ) from the National Institutes of Health, a financial support from Heidelberg Engineering, GmbH, Heidelberg, Germany, and by the LSU Biological Computation and Visualization Center. The contents presented here are solely the responsibility of the authors and do not necessarily represent the official views of the funding agencies.

The authors thank Claude F. Burgoyne, M.D., Devers Eye Institute, Portland, OR for providing us with the LEG study topographic library for evaluating our methodology.

**Madhusudhanan Balasubramanian** received his Diploma from Muthiah Polytechnic, India in 1995 and Baccalaureate of Engineering from PSG College of Technology, India in 1998, both degrees in electrical and electronics engineering. He earned his M.S. degree in systems science in 2003 and his Ph.D. degree in computer science in 2006 from Louisiana State University, Baton Rouge, LA. He is a member of IEEE since 2003.

He is currently a Postdoctoral Research Fellow at the Hamilton Glaucoma Center, University of California San Diego, La Jolla, CA. His research interests are in biomedical optical instruments for imaging microscopic structures, developing computer vision and image processing algorithms and applied mathematical and statistical approaches for studying deformation in biological structures.

Dr. Balasubramanian is also a member of ACM, ARVO, and SPIE.

**Stanislav Žabić** received his Ph.D. degree in mathematics from Louisiana State University, Baton Rouge, LA in 2005.

Dr. Žabić is currently a Staff Research Scientist at Philips Healthcare, Cleveland, OH. In 2006, he was a postdoctoral researcher in the Department of Radiology, University of Utah, investigating medical imaging techniques. His main research topics are CT image reconstruction, and mathematical modeling and simulations in physics and medicine.

**Christopher Bowd** received his Ph.D. from Washington State University in experimental psychology and neuroscience.

He is currently an Associate Research Scientist at the University of California, San Diego, Hamilton Glaucoma Center. His current work involves early detection and monitoring of glaucoma with structural imaging of the optic nerve, visual function and electrophysiological testing using standard and machine learning classifier-based analyses. He also is interested in the relationship between glaucomatous structural and functional defects, measured both cross-sectionally and longitudinally.

Dr. Bowd’s research has been funded by the National Institutes of Health, the Glaucoma Research Foundation and Pfizer Inc., and his publications include over 75 peer reviewed manuscripts and 9 book chapters.

**Hilary W. Thompson, Ph.D.** is a Professor of Biostatistics, in the School of Public Health, Louisiana State University Health Sciences Center (LSUHSC), New Orleans, LA. In addition to acting as the director of the LSU Eye Center statistical consulting and clinical trials service, he is also the statistical director of the LSUHSC Neuroscience Center biostatistics core, and the co-director of the Bioinformatics/Biocomputing core of the Louisiana INBRE grant from NIH/NCRR.

Dr. Thompson has collaborated and published extensively with computer scientists and mathematicians in research on methods for image classification and data mining methodologies.

**Peter Wolenski** received his Ph.D. in mathematics from the University of Washington in 1988.

Dr. Wolenski is currently the Russell B. Long Professor of Mathematics at Louisiana State University, Baton Rouge, LA. His research specialty is optimal control and variational analysis, and he has published 2 books and over 40 research articles.

**S. Sitharama Iyengar** is the Chairman and Roy Paul Daniels Chaired Professor of Computer Science at Louisiana State University, Baton Rouge, LA and is also the Satish Dhawan Chaired Professor at the Indian Institute of Science, Bangalore. His publications include 6 textbooks and over 380 research papers. His research interests include high-performance algorithms, data structures, sensor fusion, data mining, and intelligent systems. His research has been funded by the NSF, ONR, NASA, DoE/ORNL, U.S. Army Research Office, and the DARPA and MURI Programs.

Dr. Iyengar is a fellow of IEEE, ACM, AAAS, and SDPS. He is a recipient of IEEE awards, best paper awards, the Distinguished Alumnus award of the Indian Institute of Science, Bangalore, and other awards. He has served as the editor of several IEEE journals and is the founding editor-in-chief of the International Journal of Distributed Sensor Networks.

**Bijaya B. Karki** received his Ph.D. from the University of Edinburgh, Edinburgh, UK in 1997.

He is currently an Associate Professor in the Department of Computer Science, Louisiana State University, Baton Rouge, LA. He has been a research scholar at the University of Minnesota Supercomputing Institute and a postdoctoral researcher at the LSU Biological Computation and Visualization Center. His research interests are in the areas of scientific computing and visualization and their applications in fundamental materials research for important geophysical problems.

Dr. Karki has been awarded the prestigious CAREER grant by the National Science Foundation in 2004. He recently won the college research and university Phi Kappa Phi awards.

**Linda M. Zangwill** received her Ph.D. in epidemiology from Ben-Gurion University of the Negev, Beer-sheva, Israel.

She is currently a Professor at the Hamilton Glaucoma Center, Department of Ophthalmology, University of California, San Diego, La Jolla, CA. Her clinical research focuses on improving our ability to detect and monitor glaucomatous optic disc and retinal nerve fiber layer damage, and to understand the complex relationship between structural and functional changes in glaucoma.

Dr. Zangwill’s research has been continuously funded by the National Eye Institute of the National Institutes of Health since 1995. She has published over 150 peer reviewed manuscripts and numerous book chapters.

1. Azuara-Blanco A, Costa VP, Wilson RP. Handbook of glaucoma. Florence, KY: Taylor & Francis e-Library; 2002.

2. Kingman S. Glaucoma is second leading cause of blindness globally. Bull World Health Organ. 2004;82(11):887–888. [PubMed]

3. Quigley HA, Broman AT. The number of people with glaucoma worldwide in 2010 and 2020. Br J Ophthalmol. 2006;90(3):262–267. [PMC free article] [PubMed]

4. Lee PP, Levin LA, Walt JG, Chiang T, Katz LM, Dolgitser M, Doyle JJ, Stern LS. Cost of patients with primary open-angle glaucoma: a retrospective study of commercial insurance claims data. Ophthalmology. 2007;114(7):1241–1247. [PubMed]

5. Lee PP, Walt JG, Doyle JJ, Kotak SV, Evans SJ, Budenz DL, Chen PP, Coleman AL, Feldman RM, Jampel HD, Katz LJ, Mills RP, Myers JS, Noecker RJ, Piltz-Seymour JR, Ritch RR, Schacknow PN, Serle JB, Trick GL. A multicenter, retrospective pilot study of resource use and costs associated with severity of disease in glaucoma. Arch Ophthalmol. 2006;124(1):12–19. [PubMed]

6. Salm M, Belsky D, Sloan FA. Trends in cost of major eye diseases to medicare, 1991 to 2000. Am J Ophthalmol. 2006;142(6):976–982. [PubMed]

7. Balasubramanian M, Bowd C, Zangwill LM. Algorithms for detecting glaucomatous structural changes in the optic nerve head. In: Acharya R, Ng EY, Suri JS, editors. Image Modeling of the Human Eye. Norwood, MA: Artech House, Inc.; 2008. pp. 147–188.

8. Patterson AJ, Garway-Heath DF, Strouthidis NG, Crabb DP. A new statistical approach for quantifying change in series of retinal and optic nerve head topography images. Invest Ophthalmol Vis Sci. 2005;46(5):1659–1667. [PubMed]

9. Chauhan BC, Blanchard JW, Hamilton DC, LeBlanc RP. Technique for detecting serial topographic changes in the optic disc and peripapillary retina using scanning laser tomography. Invest Ophthalmol Vis Sci. 2000;41(3):775–782. [PubMed]

10. Sirovich L. Turbulence and the dynamics of coherent structures .1. Coherent structures. Quarterly of Applied Mathematics. 1987;45(3):561–571.

11. Belhumeur PN, Kriegman DJ. What is the set of images of an object under all possible illumination conditions? International Journal of Computer Vision. 1998;28(3):245–260.

12. Sirovich L, Kirby M. Low-dimensional procedure for the characterization of human faces. Journal of the Optical Society of America a-Optics Image Science and Vision. 1987;4(3):519–524. [PubMed]

13. Park CH, Park SN, Shin JH, Paik JK, Namkung J. Face recognition using optimized 3D information from stereo images. Image Analysis and Recognition. 2005;3656:1048–1056.

14. Gao L, Jiang J, Liang J, Wang S, Yang S, Qin Y. PCA-based approach for video scene change detection on compressed video. Electronics Letters. 2006;42(24):1389–1390.

15. Chen KW, Reiman EM, Alexander GE, Bandy D, Renaut R, Crum WR, Fox NC, Rossor MN. An automated algorithm for the computation of brain volume change from sequential MRIs using an iterative principal component analysis and its evaluation for the assessment of whole-brain atrophy rates in patients with probable alzheimer’s disease. Neuroimage. 2004;22(1):134–143. [PubMed]

16. Balasubramanian M, Iyengar SS, Beuerman RW, Reynaud J, Wolenski P. Real-time restoration of white-light confocal microscope optical sections. Journal of Electronic Imaging. 2007;16:3. [PMC free article] [PubMed]

17. Moisan Y, Bernier M, Dubois JMM. Detection of changes in a series of multitemporal ERS-1 images by principal components analysis. International Journal of Remote Sensing. 1999;20(6):1149–1167.

18. Burgoyne CF, Thompson HW, Mercante DE, Amin R. Basic issues in the sensitive and specific detection of optic nerve head surface change within longitudinal LDT TopSS images: Introduction to the LSU Experimental Glaucoma (LEG) study. In: Lemij H, Schuman J, editors. The Shape of Glaucoma. The Hague, Netherlands: Kugler Publications; 2000.

19. Nichols TE, Holmes AP. Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum Brain Mapp. 2002;15(1):1–25. [PubMed]

20. Chauhan BC, McCormick TA, Nicolela MT, LeBlanc RP. Optic disc and visual field changes in a prospective longitudinal study of patients with glaucoma: comparison of scanning laser tomography with conventional perimetry and optic disc photography. Arch Ophthalmol. 2001;119(10):1492–1499. [PubMed]

21. Artes PH, Chauhan BC. Criteria for optic disc progression with the topographic change analysis of the heidelberg retina tomograph. Invest Ophthalmol Vis Sci. 2006;47 E-Abstract 4349.

22. Bowd C, Balasubramanian M, Weinreb RN, Vizzeri G, Alencar LM, O’leary N, Sample PA, Zangwill LM. Performance of confocal scanning laser tomograph topographic change analysis (TCA) for assessing glaucomatous progression. Invest Ophthalmol Vis Sci. 2009;50(2):691–701. [PMC free article] [PubMed]

23. Artes PH, Chauhan BC. Longitudinal changes in the visual field and optic disc in glaucoma. Prog Retin Eye Res. 2005;24(3):333–354. [PubMed]

24. Wolberg G, Zokai S. Robust image registration using log-polar transform. Vancouver, BC, Canada: 2000.

25. Jones SM, Boyer AL. Investigation of an FFT-based correlation technique for verification of radiation treatment setup. Medical Physics. 1991;18(6):1116–1125. [PubMed]

26. Lewis JP. Fast template matching. Vision Interface. 1995:120–123.

27. Kirby M. Geometric data analysis: an empirical approach to dimensionality reduction and the study of patterns. New York: Wiley; 2001.

28. Dur A. On the optimality of the discrete Karhunen-Loeve expansion. Siam Journal on Control and Optimization. 1998;36(6):1937–1939.

29. Riesz F, Sz.-Nagy Béla. Functional analysis. New York: F. Ungar Pub. Co.; 1955.

30. Trefethen LN, Bau D. Numerical linear algebra. Philadelphia: SIAM; 1997.

31. Wang L, Zhang Y, Feng J. On the Euclidean distance of images. IEEE Trans Pattern Anal Mach Intell. 2005;27(8):1334–1339. [PubMed]

32. Lasko TA, Bhagwat JG, Zou KH, Ohno-Machado L. The use of receiver operating characteristic curves in biomedical informatics. J Biomed Inform. 2005;38(5):404–415. [PubMed]

33. DeLong ER, DeLong DM, Clarke-Pearson DL. Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics. 1988;44(3):837–845. [PubMed]

34. Hanley JA, Hajian-Tilaki KO. Sampling variability of non-parametric estimates of the areas under receiver operating characteristic curves: an update. Acad Radiol. 1997;4(1):49–58. [PubMed]

35. Silverman BW. Density estimation for statistics and data analysis. London; New York: Chapman and Hall; 1986.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |