PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 August 16.
Published in final edited form as:
PMCID: PMC2922066
NIHMSID: NIHMS221747

Statistics of Optical Coherence Tomography Data From Human Retina

Abstract

Optical coherence tomography (OCT) has recently become one of the primary methods for noninvasive probing of the human retina. The pseudoimage formed by OCT (the so-called B-scan) varies probabilistically across pixels due to complexities in the measurement technique. Hence, sensitive automatic procedures of diagnosis using OCT may exploit statistical analysis of the spatial distribution of reflectance. In this paper, we perform a statistical study of retinal OCT data. We find that the stretched exponential probability density function can model well the distribution of intensities in OCT pseudoimages. Moreover, we show a small, but significant correlation between neighbor pixels when measuring OCT intensities with pixels of about 5 µm. We then develop a simple joint probability model for the OCT data consistent with known retinal features. This model fits well the stretched exponential distribution of intensities and their spatial correlation. In normal retinas, fit parameters of this model are relatively constant along retinal layers, but varies across layers. However, in retinas with diabetic retinopathy, large spikes of parameter modulation interrupt the constancy within layers, exactly where pathologies are visible. We argue that these results give hope for improvement in statistical pathology-detection methods even when the disease is in its early stages.

Index Terms: Diabetic retinopathy, maximum likelihood detection, optical coherence tomography, stretched exponential distribution, visual system

I. INTRODUCTION

The retina is the first station of the brain’s visual system, and thus important for vision [1]–[4]. Therefore, almost all visual impairments, including the vast majority of types of blindness, are consequences of retinal diseases. Some of the main examples of visual diseases of the retina include diabetic retinopathy [5], [6], age-related macular degeneration [7], and retinitis pigmentosa [8]. In particular, diabetic retinopathy and macular degeneration are reaching epidemic levels [9]. Hence, all techniques that can be used in early detection of these retinal anomalies are of great interest for diagnosis and treatment of related pathologies [10]. To emphasize the importance of early diagnosis, we consider the example of diabetic retinopathy. This disease is the most common microvascular complication of diabetes [11] and one of the leading causes of blindness in the western world [12]. At the time of diagnosis, as many as 38% of type 2 diabetics manifest some degree of diabetic retinopathy [13]. Diabetic retinopathy may precede the diagnosis of diabetes by as many as seven years and progress with the duration of diabetes [14]. The appearance and progression of diabetic retinopathy are influenced by the degree of metabolic control and concomitant conditions, such as hypertension and habits such as smoking [11]. Left undiagnosed and untreated, diabetic retinopathy may progress and ultimately threaten vision. As the disease is often asymptomatic, one should detect it with periodic screening. Data from clinical trials have long shown that timely intervention with pan-retinal photocoagulation prevents vision loss due to diabetic retinopathy [15]. However, patients are not being diagnosed and treated in a timely fashion [16].

Optical coherence tomography (OCT) has recently become one of the primary methods for noninvasive mapping of the human retina [17]–[20]. Such maps can be especially useful for retinal analysis, as the retina has a strict layer organization [21]–[23]. Therefore, in each section of the map, we can see a pseudoimage of the structure of the retina. In such an image, early pathologies can often be seen as blots, i.e., local modifications of reflectance. Consequently, one may expect that automatic procedures of detection of such blots can be of help for early diagnosis of the various diseases.

The pseudoimage formed by OCT is not deterministic, but may vary probabilistically across pixels (i.e., show speckle) due to complexities in the measurement technique. There are two main sources of statistical uncertainty (noise) in OCT data. First, the medium between the OCT device and the biological tissue may contain small random variations in refraction, which tend to make the distribution of OCT intensities exponential [24]–[27]. Second, the tissue of interest is not spatially homogeneous. Different retinal cells may have different reflectances, i.e., varying indexes of refraction. Such spatial variations may destroy the exponential tendency created by the random fluctuations in the propagating medium. In the case of the retina, the random medium refractions are due to large, anisotropic molecules in anterior ocular media, such as the lens, cornea, and the vitreous humor [28]–[31]. In addition, some retinal layers contain large cell bodies of different types [21]–[23], making the retina spatially inhomogeneous. Therefore, one should not expect the distribution of OCT intensities in the retina to be exponential, but rather a mixture of exponentials. If one thinks of OCT intensities in probabilistic terms, pathological blots can be thought of as statistical perturbations of the nonpathological regions. Hence, sensitive automatic procedures of diagnosis using OCT may exploit statistical analysis of the spatial distribution of reflectance.

In this paper, we performed a statistical study of retinal OCT data. In particular, we tested various candidate probability distributions to describe the spatial variations of OCT intensity. We also analyzed interpixel dependence by studying the normalized autocorrelation function and developing a joint probability function with marginal distributions given by the stretched exponential. Furthermore, we studied the spatial dependence of the parameters of the stretched exponential both along and across retinal layers. Using the outcomes of these studies, we propose automatic statistical procedures for detection of anomalies due to retinal pathologies. We illustrate one of these procedures for a retina with diabetic retinopathy.

II. METHODS AND MODELS

A. Methods for Obtaining OCT

We used the RTVue-100 Fourier-domain OCT system (Optovue, Freemont, CA) to map the macula. A line-scan (or section) OCT image produced B-scans; we call them pseudoimages in this paper to remain consistent with terminology from computer vision from which many of our methods arise. The pseudoimages had pixels of 3 × 7.5 µm2. The finer and coarser pixel sizes were in directions perpendicular and parallel to the retinal layers, respectively. The optical resolution of the OCT parallel to the layers was about 20 µm, and thus, the system oversampled in this direction. OCT data were volumetric with 17 pseudoimages having superior–inferior retinal orientations and 17 nasal-temporal orientations. The pseudoimages corresponded to 5-mm-long sections spaced 0.25 mm apart at the central 3 × 3 mm2 and 0.5 mm elsewhere. All statistics reported in this paper come from individual pseudoimages.

As OCT data have high SNR (≈100 dB), useful visualization requires presenting intensities in a logarithm scale (see Fig. 1). Although essentially all images in the OCT literature show a logarithmic quantity, here we focus on intensity, which is the variable directly measured by the device. We will argue that the statistics of intensity provide especially useful information for the detection of retinal anomalies.

Figure 1
Examples of OCT pseudoimages of the retina and the regions used for statistical sampling. (a) Image from a nonpathological subject and sampling rectangles from the first type of method for defining homogeneous regions (Section II-D). (b) Image from another ...

The Institutional Review Board of the University of Southern California approved the OCT research protocol. In it, all subjects were treated in accordance with the Declaration of Helsinki. Informed consent was obtained from all subjects after the goals of the research and consequences of participation had been discussed. The data reported here come from 12 normal subjects and eight subjects with diabetic macular edemas (DME).

B. Statistical Model for the Distribution of OCT Intensities

As mentioned in Section I and can be seen in Fig. 1, OCT images are speckled, i.e., pixel values are distributed probabilistically. To estimate probability density functions of both OCT intensities and log-transformed values, we first constructed histograms, binning the data in 20 equally spaced intervals. For each region of interest (Section II-D), we estimated from three to ten histograms of nonoverlapping subregions. We then normalized the histograms to obtain probability density functions, and combined the data from the various subregions to estimate the mean and standard errors of the functions. The OCT intensities displayed monotonically decaying behavior.

Our goal was to find a simple model, i.e., depending on few parameters, for the OCT intensity (I) probability density functions p(I). Hence, we tried to fit several models to the data using a maximum-likelihood estimate procedure. We probed the quality of all fits with a Kolmogorov–Smirnov test (p < 0.01). The first models for p(I) that we tried were the exponential and the Weibull probability density functions. Unfortunately, these models yielded unsatisfactory fits.

A model that gave an excellent fit is the stretched exponential distribution

ps(I)=k(λ,β)e(I/λ)β
(1)

where

k(λ,β)=1λΓ((β+1)/β)
(2)

is the normalization constant of the distribution, λ > 0 and 0 < β ≤ 1 are the two parameters of the probability density function, and Γ is the gamma function.

1) Why Is the Stretched Exponential Distribution a Good Model for Retinal OCT?

Sufficiently small pixels of the pseudoimage should give rise to an exponential distribution of intensities [24]–[27]. However, because different pixels may fall in different retinal structures (e.g., different kinds of cells), the distribution in small regions will not be a pure exponential. Rather the distribution should be a mixture of exponentials of different intensity scales. A mixture of this type could give rise to an infinite number of different distributions. In our study, we attempted one class of such probability density functions, namely, the stretched exponential. We now explain the connection between mixture of exponentials and the stretched exponential distribution in more detail.

A stretched exponential distribution for a random variable x is characterized by the probability density

ps(x)={k(λ,β)e(x/λ)β,x00,x<0
(3)

where the normalization factor k(λ,β) is as in (2). For β = 1, (3) reduces to the exponential density, but for β < 1, the equation shows an elongated (stretched) tail. If one calculates the moment of order k of (3), the result is

E(xk)=Γ(k+1/β)Γ(1/β)λk,  k=1,2,.
(4)

Hence, in particular, the mean, second moment, variance, and coefficient of variation are

E(x)=Γ(2/β)Γ(1/β)λ
(5)

E(x2)=Γ(3/β)Γ(1/β)λ2
(6)

σ2(x)=Γ(3/β)Γ(1/β)Γ2(2/β)Γ2(1/β)λ2
(7)

and

CV=(Γ(3/β)Γ(1/β)Γ2(2/β))1/2Γ(2/β)
(8)

respectively. Inspection of (5) and (6) shows a one-to-one correspondence between pair (E(x), E(x2)) and (λ,β). In turn, (8) evidences that the coefficient of variation is independent of λ.

To understand the connection between the stretched exponential distribution and the mixture of exponentials, we begin with the notion of failure rate. For a distribution with support in [0,∞) and with probability density p(x), the failure rate is defined as

f(x)=p(x)xp(ξ)dξ.
(9)

In particular, for a stretched exponential distribution, the failure rate is

fs(x)=e(x/λ)βxe(ξ/λ)βdξ
(10)

which is nonincreasing with respect to x [set membership] [0,∞). Therefore, this distribution belongs to the class of the decreasing-failurerate (DFR) distributions [32]. In turn, DFR distributions are intimately connected to the class of mixture-of-exponentials distributions [33]. The following theoretical results hold for these distributions [32], [33].

Theorem 1

A mixture of exponential distributions, i.e., a distribution characterized by the density

p(x)=01τe(x/τ)µ(dτ)
(11)

where µ is a probability measure on [0,∞), is DFR.

Moreover, on the basis of results proved by Langberg et al. [33], the following statement holds.

Theorem 2

A DFR distribution with a continuously differentiable failure rate is a mixture of exponential distributions.

From the earlier results, it follows that a stretched exponential distribution is a mixture of exponentials

ps(τ)=01τe(x/τ)µs(dτ)
(12)

with the measure µs determined by the pair (λ,β). Equation (12) allows connecting moments of a stretched exponential distribution to the corresponding moments of µs as follows:

E(xk)=00xkτe(x/τ)dxµs(dτ).
(13)

Hence, using (4), one gets

Γ((k+1)/β)Γ(1/β)λkk!=0τkµs(dτ).
(14)

Berberan–Santos et al. [34] performed a general investigation about the mixture-of-exponentials representation of the stretched exponential distribution. They stressed that µs(dτ)/τ defined in (12) plays the role of the inverse Laplace transform of ps. Moreover, the numerical computation of µs given (λ,β) was discussed along with the related difficulties.

C. Estimation of Autocorrelation

To study the statistical dependence between nearby pixels in pseudoimages, we computed the sampled spatial normalized autocorrelation function. This was done by first estimating the autocorrelation function. We estimated this function in square windows of 21 × 21 pixels, corresponding to maximal distances of ≈ 170 µm. Then, we normalized it by the variance, i.e., our normalized autocorrelation values were numbers from −1 to 1. We also tested the probability that each normalized autocorrelation value was different from zero using the Fisher’s z statistic [35] (p < 0.001).

1) Model for the Autocorrelation

This section focuses on the mathematical aspects of the model for autocorrelation, leaving physical interpretations to the next Section II-C2 that the stretched exponential distribution gave a good fit to the OCT-intensity data and constrained our search of a model for the autocorrelation function. We had to look for a 2-D distribution of (possibly) nonindependent variables having stretched exponentials as 1-D marginals. The model that we came up with was based on a mixture of a product of exponential distributions. The representation as mixture of exponentials (12) provided a way to find such a 2-D joint distribution. We considered the mixture of a product of exponential distributions, i.e.,

p(x,y)=001τ11τ2e(x/τ1)e(y/τ2)µ(dτ1,dτ2).
(15)

This mixture depends on a weight function µ, which one may understand as follows. As discussed in Section II-B, although individual pixels should yield exponential distributions, a spatially distributed set of pixels may not, since, for example, different pixels may come from different cellular processes. The weight function measures the probability density for pixels being in the same cellular process. In other words, the development of the weight function is based on a local space-homogeneity assumption.

We developed the weight function µ to ensure the stretched exponential marginals. The density in (15) has marginals

px(x)=01τ1e(x/τ1)µx(dτ1)
(16)

and

py(y)=01τ2e(y/τ2)µy(dτ2)
(17)

where µx and µy are the marginals of µ. Thus,px and py will be of the stretched exponential type as soon as µ is such that its marginals, µx and µy, satisfy (12), with the corresponding values of (λ,β), namely (λx, βx) and (λy, βy).

We now develop a procedure to obtain such an appropriate 2-D distribution µ. Let us introduce three numbers [var phi], ρ, and ω, such that

0φ,ρ,ω1.
(18)

Then, given two pairs of nonnegative numbers λ1 ≠ λ2 and β1 ≠ β2, we define the pairs (λxx) and (λyy) as follows:

λx=(2ω)λ1+ωλ22,βx=(2ω)β1+ωβ22λy=ωλ1+(2ω)λ22,βy=ωβ1+(2ω)β22.
(19)

Hence, ω controls how much µx and µy define the same stretched exponential distribution. If ω = 0, then βx = β1, βy = β2, λx = λ1, and λy = λ2, making the distributions different. As ω increases, βx and βy, and λx and λy start moving toward each other. In the limit of ω = 1, βx = βy and λx = λy, yielding the same distribution. Next, we choose µ in (15) as the following measure on [0, + ∞)2:

µ(dτ1,dτ2)=ρ2[µx(dτ1)+µy(dτ2)]δ(τ1τ2)+µx(dτ1)µy(dτ2)ρ2[µx(dτ1)µy(dτ2)+µx(dτ2)µy(dτ1)]
(20)

where δ denotes the Dirac delta function. One can show that µ is a probability measure provided that

ρ21+Axµy(dτ1)Ayµx(dτ2)Axµx(dτ1)Ayµy(dτ2)
(21)

for any Borelian sets Ax and Ay in [0,∞), and one can show that the marginals of µ are µx and µy. The parameter ρ controls how correlated the variables of µ are. If ρ = 0, then from (20),µ(τ12) = µx1y2), making τ1 and τ2 independent. The larger ρ gets, the more the δ function dominates (20), making the occurrence of τ1 = τ2 more likely. Therefore, 1 – ρ is the probability that we pick τ1 and τ2 independently.

Besides defining µ with a correlation structure, we also allow correlation to occur at the level of the measurements x and y. For this purpose, we use the µ given by (20) to define the probability density function [p with macron](x,y) as in (15), with marginals px(x) and py(y). From these functions, we define the joint density

p(x,y)=φ2[px(x)+py(y)]δ(xy)+p¯(x,y)φ2[p¯(x,y)+p¯(y,x)].
(22)

As with µ, one can show that p(x,y) is a probability density with marginals px(x) and py(y), provided that

φ21+Axpy(x)dxAypx(y)dyAxpx(x)dxAypy(y)dy.
(23)

The parameter [var phi] controls the correlation at the level of the measurements themselves. If [var phi] = 0, then from (22) and the definition of [p with macron](x,y), p(x,y) arises directly from µ through (15). The larger [var phi] gets, the more the δ function dominates (22), making the occurrence of x = y more likely. Consequently, we can interpret 1 – [var phi] as the probability that we pick x and y independently.

We next proceed to calculate the covariance K(x,y) of this model. By definition, this function is

K(xy)=E(xy)E(x)E(y).
(24)

To calculate E(xy), we start from (22) to write

E(xy)=00ξηp(ξ,η)dξdη.
(25)

We then expand [p with macron], px, and py by using (15)(17) with µ, µx, and µy given by (20). Finally, we expand the results in terms of Γ functions by using (5), (6), and (14). The result is

E(xy)=(φ2+(1φ)ρ4)[λx2Γ(3βx)Γ(1βx)+λy2Γ(3βy)Γ(1βy)]+(1φ)(1ρ)λxλyΓ(2βx)Γ(2βy)Γ(1βx)Γ(1βy).
(26)

Substituting this equation for E(xy) in (24), and using (5) for E(x) and E(y), we obtain for the covariance

K(x,y)=(φ2+(1φ)ρ4)[λx2Γ(3βx)Γ(1βx)+λy2Γ(3βy)Γ(1βy)](φ+ρ(1φ))λxλyΓ(2βx)Γ(2βy)Γ(1βx)Γ(1βy).
(27)

Dividing this equation by the product of the two standard deviations [square roots of (7)] gives the correlation coefficient shown in (28) at the bottom of this page.

R(x,y)=(φ2+(1φ)ρ4)[λx2Γ(3βx)Γ(1βy)+λy2Γ(3βy)Γ(1βx)](φ+ρ(1φ))λxλyΓ(2βx)Γ(2βy)λxλy(Γ(3βx)Γ(1βx)Γ2(2βx))(Γ(3βy)Γ(1βy)Γ2(2βy))
(28)

In this paper, we attempted to make our measurements in homogeneous zones of the retina, i.e., with ω = 1. In the rest of this paper, we thus focus on this condition. If ω = 1, from (19), we have

λx=λy=λ=λ1+λ22  βx=βy=β=β1+β22
(29)

and thus, µx = µy and px = py, forcing (21) and (23) to yield ρ,[var phi] ≤ 1. In addition, (21) simplifies to

R(x,y)=φ+(1φ)ρ2Γ(3β)Γ(1β)2Γ2(2β)Γ(3β)Γ(1β)Γ2(2β).
(30)

This correlation coefficient has five noteworthy properties. First, it does not depend on λ, thus resembling the behavior of the coefficient of variation [see (8)]. Second, because Γ(3/β)Γ(1/β) ≥ 2Γ2(2/β), R(x,y) ≥ 0. This should not be surprising because (20) and (22) impose only positive correlations by creating instances where τ1 = τ2 and x = y. Third

limφ1R(x,y)=1.
(31)

i.e., as expected, if the probability that the two measurements are equal goes to 1, so does the correlation. Fourth, if [var phi] → 0, then R becomes proportional to ρ, and if ρ → 0, then R[var phi]. In other words, in the absence of one source of correlation [imposed by (20) or (22)], the other source dominates. A corollary of this conclusion is

limφ0ρ0R(x,y)=0
(32)

i.e., if one eliminates both sources of correlation in the model, the correlation goes to zero. Fifth, if β = 1, i.e., the marginals are exponential, then R(x,y) = [var phi], because Γ(3)Γ(1) = 2Γ2(2) Thus, in particular

limφ0β1R(x,y)=0.
(33)

2) Interpretation and Specification of the Model

A physical interpretation of the model in the context of our study is as follows. Variables x and y represent OCT intensities. Their individual distributions are stretched exponentials [(3) and (4)], and these intensities may exhibit correlation [(28) or (30)]. Typically, the two intensities are measured at two different pixels separated by a given distance vector r. We expect that, as the magnitude of r increases, the value of the normalized autocorrelation function should decrease, being 1 when the distance is zero and converging to 0 when the distance is very large. Because parameters ρ and [var phi] control correlation, [(28) or (30)], we propose that they fall with |r|. For example, we can rewrite (30) as

R(r)=φ(r)+(1φ(r))ρ(r)2Γ(3β)Γ(1β)2Γ2(2β)Γ(3β)Γ(1β)Γ2(2β).
(34)

This fall is consistent with the interpretation of the parameters given after (21) and (23). In particular, ρ should decrease with |r|, as the measurement would be made from different cellular processes. Similarly, [var phi] should decrease with |r|, since the medium between the OCT device and the biological tissue contains spatially independent random variations in refraction. Furthermore, the image would have to be still and perfectly focused to prevent loss of independence between neighbor intensities. The loss of independence by joint refraction, eye motion, or loss of focus should diminish with distance, causing [var phi] to fall with |r|. According to these physical interpretations for [var phi] and ρ, we propose that the former falls faster with |r| than the latter, because focusing and eye-motion errors are minimal and refraction fluctuations would be in molecular spatial scales, which are smaller than cellular ones. For simplicity, we model the spatial dependences of [var phi] and ρ as the exponential decays

ρ(r)=e|r|/υρ  φ(r)=e|r|/υφ
(35)

where υρ > υ[var phi] are length constants. These dependences ensure that the normalized autocorrelation function converges to 1 at short distances [(31)] and to 0 at large distances [(32)].

We can also give physical interpretations to the other three parameters of the model, namely, ω,λ, and β. The number ω controls the distance between the marginals px and py, which coincide for ω = 1. In other words, for ω < 1, the two intensities come from cell populations with different distributions. For instance, this happens whenever the two pixels belong to two statistically different layers [e.g., outer-nuclear layer (ONL), and inner segment (IS)/retinal pigment epithelium (RPE)]. Our efforts to work with homogeous layers thus justify our decision to work with ω = 1 [(30) instead of (28)].

The parameter λ is an intensity scale [(8)], with bearing on neither the shape of the distribution [(8) and (29)] nor the correlation structure of the model [(30)] when ω = 1. In contrast, β controls both the shape of the distribution and its correlation structure. An especially interesting situation occurs for exponential marginals, i.e., for β = 1. If [var phi] is small (which is the typical case for us), then we obtain R(x,y) = 0 [(33)]. This equality will lend support to our model, as we will show that outside the retina and in the ONL, R(x,y) = 0 for |r| > 0. We will also show that the OCT intensity distribution in these zones is nearly exponential, i.e., β ≈ 1. From the perspective of our model, stretched exponential distributions with β = 1 are interesting because for them, the function µ is a delta function [(12)]. This means that all the measurements come from the same exponential distribution. In other words, the retinal region where we are making the measurements does not contain intensity fluctuations due to cellular processes. In conclusion, the model suggests that for regions with few cells (e.g., humor vitreous), intensities are exponentially distributed and without spatial correlations. The same should hold for regions with densely packed cells, which are much smaller than measurement pixels (e.g., the ONL).

D. Obtaining Regions for Statistical Analysis

Our goal was to construct regions in which the probability distribution of OCT intensities varied as little as possible. In this study, we tried two types of methods for defining such regions. The first attempted to find homogeneous regions by sampling small pieces of the retina. The second took advantage of the histological homogeneity along retinal layers. The retina has a multilayer spatial organization. Some of these layers are as follows: nerve–fiber layer (NFL), ganglion-cell layer (GCL), inner plexiform layer (IPL), inner nuclear layer (INL), outer plexiform layer (OPL), outer-nuclear layer (ONL), inner segments of photoreceptors (IS), and retinal pigment epithelium (RPE). In more detail, the two types of methods are as follows.

  1. We used small rectangles to sample regions that appeared similar by eye [see Fig. 1(a)]. Rectangles were parallel to the retinal layers, with sizes ranging from 20 × 10 pixels(≈ 150 × 30 µm2) to 50 × 30 pixels(≈ 375 × 90 µm2).
  2. We developed algorithms for detecting four or five major borders in the image. They included the retinal-vitreous border, the border between the NFL and the GCL/IPL/INL/OPL regions, the border between the GCL/IPL/INL/OPL and the ONL, the border between the ONL and the IS/RPE regions, and the lower border of the RPE with the choroid see Fig. 1(b)] [Fig. 1(d) presents our nomenclature for the various regions]. One can use OCT to subdivide the normal retina into more than the four major regions (NFL, GCL/IPL/INL/OPL, ONL, and IS/RPE) that we used here [17], [36]. However, we preferred to use these four regions, as certain disorders, such as DME, fuse some of the layers in terms of OCT imaging [37]. By keeping a coarse subdivision, we could compare DME and normal retinas fairly. Any way, our statistical study below used many more than the four major divisions shown in Fig. 1. We used up to 20 statistically homogeneous divisions. With such large number of layers, we ensured that our conclusions were not limited by a poor subdivision of the retina.

To find these four or five borders, all algorithms broadly measured zones of statistically constant reflectance and then detected sharp gradients between these zones.

E. Methods for Estimating Distribution Parameters and Their Spatial Variations

In this study, we performed fittings of (1), without constraining λ and β. We observed that their optimal values were always positive, and although the best β was sometimes larger than 1, it was not statistically significantly so. As shown by (5) and (6), there is a one-to-one correspondence between λ and β, and the first two moments of the distribution in (1). Consequently, one can obtain a first estimate of λ and β directly from estimates of the mean and the variance. We use this first estimate of these parameters as initial condition in our maximum-likelihood estimate of λ and β. In this study, we estimated these parameters from samples obtained in small regions as defined in Section II-D.

We hypothesized that λ and β would be substantially constant in narrow regions parallel to the retinal layers. In contrast, λ and β should vary strongly in perpendicular orientations. To test this hypothesis, we used the two types of methods described in Section II-D. With type 1, we placed from three to seven rectangles as far as possible apart inside individual retinal layers identified by eye [see Fig. 1(a)]. With type 2, we divided the small closed regions [see Fig. 1(b) and Fig. 1(c)] in 75-µm-wide nonoverlapping samples. We then estimated their median OCT intensities. To the spatial dependence of these medians, we fitted third- to fifth-order polynomials to with the random sample consensus (RANSAC) method [38]. The intensity in each pixel was then multiplied by the ratio between the polynomial value at the center of the sample and its median. This procedure removed small inhomogeneities across space due to imaging artifacts. Next, we estimated λ and β [(1)], as well as k [(2)] for each sample, probing how they varied along each small, closed region [see Fig. 1(b) and Fig. 1(c)]. We found that the estimation of λ and β was somewhat ill conditioned, because their estimates exhibited some positive correlation even along homogeneous retinal layers. We thus adopted a procedure that regularized the estimation in the direction parallel to the retinal layers. The regularized procedure had four steps, obtaining: 1) unconstrained fits over a layer; 2) the median β (β*) of those fits; 3) the spatial dependence of λ (λ*(x)) with the constraint β = β*; and 4) the spatial dependence β(β*(x)) with the constraint λ = λ*(x).

III. RESULTS

A. Distribution of OCT Intensities

In this section, we focus on the probability density functions for the OCT intensity. An example for a normal subject can be seen in Fig. 2; we obtained similar distributions for 12 nonpathological and eight DME subjects. Fig. 2 shows that the distributions of OCT intensities inside the sampling rectangles are statistically similar when they are contiguous and inside a single retinal layer. The distributions are monotonically decaying. In other words, the most likely OCT intensities are low, whereas high intensities are rare.

Figure 2
Probability density functions of OCT intensities. The data for this example come from the same image as in Fig. 1(a). However, for this figure, we used three sampling rectangles (375 × 60 µm2). They were contiguous and nonoverlapping in ...

As explained in Section I and Section II-B, the exponential distribution would be the simplest, physically justifiable model for the decaying distribution of OCT intensities. An interesting alternative is the Weibull distribution. This distribution is often used in the field of life data analysis due to its flexibility [39]. Unfortunately, neither the exponential nor the Weibull fit was completely satisfactory. These models failed, because the data distribution had a long tail. We confirmed the statistical inadequacy of these models with a Kolmogorov–Smirnov test (p < 0.01). The small but significant failures of these models were true for all our subjects (N = 20). (However, as we will discuss in Section IV-B, the ONL lends itself to a good exponential fit. In addition, one can obtain good exponential fits for zones outside both the retina and the RPE—Section IV-A.)

Because neither the exponential nor the Weibull distribution provides good fits over the entire retina, we tested the applicability of the stretched exponential probability density function [(1)], following the rationale provided in Section II-B. There we reasoned as follows: because different pixels may fall in different retinal structures (e.g., different kinds of cells), the distribution in small regions will not be a pure exponential. Rather the distribution should be a mixture of exponentials of different intensity scales (Section II-B1). One such distribution much used in Physics is the stretched exponential. We thus attempted fitting stretched exponentials to the data (see Fig. 3).

Figure 3
Fits of the data with the stretched exponential model. (a) Fit of the data of Fig. 2 [(1);λ = 1200; β = 0.869]. (b) Fit in a nonpathological ONL, using a new subject (λ = 302; β = 1.16−β was not statistically ...

Fig. 3 shows that the stretched exponential probability density function provides an excellent fit to the data. The good quality of the fits holds whether we are in a normal [see Fig. 3(a)–(c)] or a DME patient blot [see Fig. 3(d)] zone of the cellular tissue. (The larger error bars shown in Fig. 3(d) are due to the DME blots being relatively small, and thus yielding small amounts of data.) Furthermore, fits are good regardless of the zone of the OCT image where one makes the measurement. Consequently, fits captured fast decays (small λ; e.g., ONL and DME patient blots), slow decays (large λ; e.g., NFL and IS/RPE), long tails (small β; e.g., IS/RPE), and short tails (large β; e.g., ONL). Using a Kolmogorov–Smirnov test (p < 0.01), we determined that all the fits shown in Fig. 3 were statistically undistinguishable from the data. The same held for 97% of the data that we probed in all our subjects and in different zones of OCT images. For the other 3%, the failure of the fit was marginal, i.e., with p just barely above 0.01.

In conclusion, the stretched exponential probability density function can model well the distribution of intensities in OCT pseudoimages.

B. Autocorrelation in OCT pseudoImages

The stretch exponential describes the distribution of OCT intensities, but says nothing about the dependency across neighbor pixels in a pseudoimage. To quantify this dependency, we computed the spatial normalized autocorrelation of OCT intensities in small rectangles (Section II-C). Examples for two retinas are shown in Fig. 4.

Figure 4
Normalized autocorrelation functions for two small retinal regions. The spatial profile appears in a false color scheme with 1 being in red, with 0 being in dark blue, and in between correlations appearing in brighter tones of blue. (a) This autocorrelation ...

The results shown in Fig. 4 show two kinds of small, but significant correlation between neighbor pixels. The correlation in the horizontal direction (i.e., parallel to the retinal layers) has a long range. Such a long-range correlation in this direction is not surprising. The optical resolution of the system in this direction is 20 µm (Section II-A). Much more interesting is the shorter range correlation in the vertical direction (i.e., perpendicular to the retinal layers—A-scans). In this direction, we found relatively strong correlations only when measuring OCT intensities with pixels of around 5 µm [see Fig. 4(a) and (b)]. Cross sections around the center of the normalized autocorrelation function reveal that within 5 µm the correlation is typically around 0.3 [see Fig. 4(c) and (d)]. This short-distance correlation holds for all retinal layers except the ONL, which only shows correlation at zero distance. (A discussion of the ONL will appear in Section IV-B.) Over larger distances, the behavior of the autocorrelation varies over different retinal regions. We found that ≈90% of all regions displaying nonindependent neighbor pixels show a fast spatial decay after about 5 µm [see Fig. 4(a) and 4(c)]. For the other ≈10%, a slowed decay was observed after the initial 5 µm [see Fig. 4(b) and (d)]. For example, in Fig. 4(b) and (d), one observes statistically significant positive correlations for distances of up to 12 µm.

Section II-C1 discusses two possible reasons for the observed correlation. We then use them to develop a correlation model compatible with the stretched exponential distribution of intensities. The behavior of the normalized autocorrelation in this model is described by (34) and (35). The model identifies two sources of correlation. They are captured by the functions ρ and [var phi], and their associated parameters υρ and υ[var phi]. The function ρ expresses correlation due to different regions of the retina having different mean OCT intensities because of reflectance from different cellular processes. In turn, the function [var phi] expresses correlation due to eye motion or loss of image focus, or locally joint refraction in the medium between the OCT device and the biological tissue. This latter three mechanisms are expected to cause correlation that falls more rapidly with distance than the earlier one. Because υρ and υ[var phi] are the lengths constant in our model, this expectation imposes υρ > υ[var phi]. As Fig. 4(c) and (d) shows, fitting the data by optimizing υ[var phi] and constraining υρ[var phi] = 5 yields excellent results. Moreover, the optimal fitting parameters are consistent with known values of ganglion-cell soma sizes, which range from 5 to 25 µm in mammalian retinas [40]–[43]. For all our subjects (N = 20 subjects), the fit was good and the range of optimal υρ was between 10 and 27 µm.

In contrast, fitting the data with lower values of υρ[var phi] yielded an underestimation of correlation at distances of 6–12 µm in data as that of Fig. 4(d). A similar underestimation occurred when replacing the exponentials in (35) with Gaussians. Similarly, using just one of the two sources of correlation for the fit did not yield good fits. Zeroing ρ or [var phi] yielded underestimation of correlation. For instance, for Fig. 4(c), when [var phi] = 0, the largest normalized correlation the model predicts for two separate pixels is 0.098 (β = 0.815; ρ → 1). Yet, we measured correlations of 0.32 and 0.25 at the neighbor directions perpendicular and parallel to the retinal layers, respectively. Similarly, we predicted a maximal two-pixel correlation of 0.068 for Fig. 4(d), but observed 0.35 and 0.33 in the perpendicular and parallel directions, respectively. Although [var phi] alone could yield better estimates than ρ alone, using only [var phi] still yielded underestimated correlations at distances between 9 and 12 µm.

Hence, one cannot assume independence in the analysis of intensities in neighbor pixels. However, a simple joint probability model (Section II-C-1) and its associated autocorrelation may be used to capture the observed dependence.

C. Space Dependence of OCT Intensity Distribution

Fig. 3 showed that the stretched exponential distribution provides an excellent fit to the distribution of OCT intensities. How do the parameters of the fits vary with position in directions parallel and perpendicular to the retinal layers? To answer this question, we studied the spatial dependence of these parameters both along and across retinal layers (Sections II-D and II-E). We began by looking at the variation of median OCT intensity and the covariation between estimated λ and β [(1)]. Results on the median appear in Fig. 5.

Figure 5
Spatial dependence of the median OCT intensity. These data were obtained from the same retina as in Fig. 3(c). Median OCT intensity appears (in semilogarithmic axes) as function of retinal position along (lower abcissas) and across (higher abscissas) ...

Fig. 5 shows that the median OCT intensity is relatively constant along the RPE of a normal retina. The same was observed over for all retinal layers in all normal retinas studied. In contrast, the median OCT intensity varied widely across retinal layers. The relative difference between the layer with lowest median (ONL) and that of highest median (RPE) was almost two orders of magnitude. The median was medium at the NFL (i.e., close to the optic-fiber layer), fell throughout the GCL/IPL/INL/OPL, reaching a minimum at the ONL (i.e., at ≈ 200 µm from the optic-fiber layer). Then, the median rose sharply at the IS/RPE. The median patterns shown in Fig. 5 hold for all retinas studied.

That the relative changes in median OCT intensity along individual layers were small gave us hope. Perhaps the parameters of the stretched exponential fits would also be relatively constant along layers. We thus studied the dependence of the parameters of the stretched exponential fit in the direction parallel to the retinal layers [see Fig. 6(a)–(c)]. In all layers, λ, β, and k (2) were relatively constant. For example, in Fig. 6, the coefficients of variation for λ, β, and k were only 0.19, 0.014, and 0.22, respectively, in the IS/RPE. Respective coefficients were 0.22, 0.065, and 0.21 in the ONL, 0.14, 0.052, 0.12 in the GCL/IPL/INL/ONL, and 0.38, 0.083, and 0.42 in the NFL. (We observed similar values in 11 other retinas.) In contrast, these parameters changed substantially across layers [see Fig. 6(d)–(e)]. The parameter λ was lowest at the ONL (λ = 170 ± 30; mean ± standard error over 12 normal retinas). This parameter then increased at the GCL/IPL/INL/OPL (λ = 300 ± 60). Finally, λ reached its highest values at the extreme layers, i.e., at the NFL (λ = 940 ± 140) and, specially, at the IS/RPE (λ = 2490 ± 400). In contrast, β varied much less. Its minimum and maximum values were at the IS/RPE (0.69 ± 0.07) and ONL (1.07 ± 0.21), respectively. The NFL (0.81 ± 0.12) and the GCL/IPL/INL/OPL (0.78 ± 0.08) yielded intermediate values. Because k was a parameter dependent on λ and β [(2)], and because λ varied much more than β, k had almost, but not precisely, an inverse relationship with λ [compare Fig. 6(d) and (f)].

Figure 6
Dependence of the parameters of the stretched exponential fit in directions both parallel (a)–(c) and perpendicular (d)–(f) to the retinal layers. The retina was the same as in Fig. 3(c). In (a)–(c), fits were to data taken over ...

Because fit parameters are relatively constant along retinal layers, we conclude that the stretched exponential distribution may be useful to detect anomalies that may affect localized regions of a layer.

D. OCT Intensity Distribution in Diabetic Retinopathy

To explore a little deeper the possibility of using the stretched exponential distribution to detect retinal anomalies, we began measuring its parameters in patients with DME. We focused on blots of DME retinas, as by eye, they appeared significantly deeper than in normal retinas. Blots are not necessarily the main anomaly in DME retinas, which can also present fusion and thickening of layers [37]. However, blots are specially amenable to the techniques we developed here. (We plan to expand efforts [44] to study the statistics of retinal-layer thickening in a future study.)

We already knew that the stretched exponential distribution could fit well data in blots of patients with DME [see Fig. 3(d)]. However, we were not certain whether the modulation of parameters due to the disease would be sufficiently deep for detection. After all, parameters are somewhat variable even in normal individual layers [see Fig. 6(a)–(c)]. We thus repeated the parametric estimation of Fig. 6(a)–(c) in patients with DME. Results for the same patient as in Fig. 1(c) and (d) appear in Fig. 7.

Figure 7
Dependence of the parameters of the stretched exponential fit as a function of position in the RPE of a patient with diabetic retinopathy. The retina was the same as in Fig. 1(c) and (d). Fits were to data taken over the entire breadth of the RPE, with ...

The results of Fig. 7 show that DME causes spikes in both the λ and k dependences on retinal position. In contrast, we do not observe clear signatures of the edemas with the β parameter [see Fig. 7(b)]. Fig. 7(a) and (c) has four spikes, each coinciding with the edemas apparent in Fig. 1(c). These spikes in the λ curve are downward [see Fig. 7(a)]. This is not surprising given that blots in DME patients cause a reduction of OCT intensities ([see Fig. 1(d)]. Because the mean of the OCT intensity distribution increases with λ [(5)], it thus should fall in the edemas. In turn, the edema spikes in the k curve were upward [see Fig. 7(c)]. As before, the reason for the k spikes was that it depended on λ and β [(2)], and that λ varied much more than β. Consequently, k should have had an almost, but not precise, inverse relationship with λ [(2)]. Interestingly, however, the variability of the baseline of the k curve [see Fig. 7(c)] was much smaller than that of the λ curve [Fig. 7(a)]. In addition, the k curve showed great sensitivity to the anomalies of the DME retina. For instance, the coefficients of variation of normal k curves were only around 0.2 in the RPE (Section III-C). In contrast, the four spikes shown in Fig. 7(c) were 2.7, 1.8, 11.2, and 3.5 times larger than the baseline. Therefore, these spikes corresponded to SNRs of about 9, 4, 50, and 13, respectively. Such large SNRs give us hope that this parameter may provide a good tool for detection of anomalies in diabetic retinopathy and perhaps other retinal diseases.

IV. DISCUSSION

A. What Does a Stretched Exponential Distribution of OCT Intensities Mean?

One of the main findings in this paper is that, among many candidate distributions, the stretched exponential probability density function accounts especially well for the retinal OCT intensity distribution (see Fig. 3). This function has been used before successfully to model other biological processes [34], [45], [46]. Why does this probability density function fit our data so well? To answer this question, we began with the result in Section II-B1. It showed that the stretched exponential is an example of a mixture of exponential distributions. To explain how this mixture arises, we proposed a physical model, which were compatible with our knowledge of both OCT and the retina.

We already addressed our physical model for the stretched exponential in Section I. In that model, small random variations in refraction in the medium between the OCT device and the biological tissue causes the tendency of OCT intensities to exhibit exponential distributions. These variations are due to large, anisotropic molecules in anterior ocular media [28]–[31]. They cause the phase of individual rays of light to undergo a random walk [24]–[27]. Consequently, because the number of random steps is large, the final amplitudes of coherent light have a Gaussian distribution (due to the central limit theorem). This Gaussian distribution of amplitudes becomes an exponential distribution of OCT intensities, since the intensity is the square of the amplitude. However, the distribution of retinal OCT intensities is not exponential, but a stretched exponential. Why is the exponential distribution stretched? In Section II-B, we suggested that the reason was that different cells in the retina and the pigment epithelium have different reflectance, yielding a distribution of OCT intensities, which is a mixture of exponentials (Section II-B1).

We call this the cellular-origin hypothesis for the retinal OCT stretched exponential distribution of intensities. This paper provides no direct test for this cellular-origin hypothesis. However, several lines of data are consistent with it. For instance, if it were true, one would expect that outside the cellular tissue, the distribution of OCT intensity should be purely exponential. To test this idea, we sampled these intensities in the vitreous humor. Fig. 8(a) shows that the exponential distribution fits the vitreous-humor OCT intensity data well. As before, we confirmed that this distribution is not statistically different from that of the data with a Kolmogorov–Smirnov test (p < 0.01). And as before, the goodness of the fit of this distribution applies to all our subjects (N = 20). Therefore, this observation is consistent with the cellular-origin hypothesis for the stretched exponential distribution.

Figure 8
Distribution of OCT intensities in the vitreous humor and evidence of cellular processes in OCT images. The data in this figure come from the same retina as in Fig. 3(c). (a) Distribution of intensities in the vitreous humor obtained with the same methods ...

B. What Does Nonzero Autocorrelation Between OCT Intensities Mean?

Another important finding is that for the RPE and most retinal layers, OCT images display nonzero autocorrelation for distances between 3 and over 20 µm. We developed a model that accounted for these autocorrelation data (Section II-C). Appropriately, this model had a stretched exponential marginal distribution. The model had two sources of correlation. The first was based on the cellular hypothesis discussed in Section IV-A. In turn, the second source of correlation could be related to other physical factors such as eye motion or loss of image focus, and as locally joint refraction in the medium between the biological tissue and the OCT device (Section II-C).

An interesting prediction of the cellular-origin hypothesis is that cellular structures should contribute to spatial autocorrelation in the retina and the pigment epithelium, but not outside them. The reason for this contribution to autocorrelation is that our sampling pixels (3 × 7.5 µm2) are smaller than many cells in the retina [40], [41], [43]. Therefore, different local sets of measurement may lead to local correlation (see Fig. 4). Fig. 8(b) confirms the prediction that in the absence of cells, we find little or no correlation. Further confirmation comes from similar results in the ONL (Section III-B). In this layer, somas are small [47]–[49] (mean diameter in central retina < 2.5 µm), and thus, according to our hypothesis, the ONL should both exhibit little or no autocorrelation, and display an exponential distribution of intensities [see Fig. 3(b)]. Equation (33) confirms the link between exponentiality and zero autocorrelation, because if β = 1 (and [var phi] = 0), then R = 0 (see also Section II-C2). (The reason for the necessity of the condition [var phi] =0 is that [var phi] ≠ 0 indicates the presence of the second source of correlation mentioned before.)

One gets further support for the cellular-hypothesis explanation of correlation if one magnifies images such as those shown in Fig. 1. For example, as shown in Fig. 8(c), one does not see spatial structures in the ONL (under the arrows). However, in the INL (yellow regions indicated by arrows), one can appreciate regions where if a pixel has higher intensity than average, then neighbor pixels also tend to be more intense. Therefore, these regions, which are not apparent in the ONL contribute to nonzero autocorrelation in the INL. Interestingly, the sizes of these regions shown in Fig. 8(c) are around 20–30 µm. These sizes coincide with some of the largest soma sizes in the retina [40]–[43]. Consequently, the B-scan may occasionally cut a soma at its largest diameter, yielding structures as in Fig. 8(c). These diameters also coincide with the largest values of υρ in the fits of Fig. 4. This parameter is important because it represents the spatial scale of the cellular source of correlation [(34) and (35)].

C. Existence of Statistically Homogeneous Layers

Another important finding in this paper is that the retina and the pigment epithelium show a layered structure in terms of the statistical parameters of the OCT intensity distribution. These parameters are substantially constant along these layers [see Fig. 5 and Fig. 6(a)–(c)], but not across them [see Fig. 6(d)–(f)]. Consequently, the OCT image conserves the layered structure found by histological studies of the retina [50], [53], [54]. Other studies proposed a cellular identification of OCT layers [17], [36], [44], which we use throughout this paper (see Fig. 1). The existence of such layers greatly increases our hopes of finding tools for early diagnoses of retinal pathologies, as we will address in Section IV-D.

A quantitative analysis within and across OCT-defined retinal layers (Section II-D) reveal some interesting behaviors. In all retinal OCT images, the IS/RPE and the NFL are always the highest reflecting layers, whereas the ONL always reflects the least. The fitted β of the stretched exponential is relatively constant [see Fig. 6(e)]. Therefore, the modulation of reflectance across layers causes a similar modulation in the fitted λ [see Fig. 6(d)]. This parameter varies by about 1.5 log units across retinal layers. Although β varies much less, it shows interesting modulation. As aforementioned, β peaks at the ONL, reaching 1 there. The lowest values of β are in the IS/RPE, where they can fall to near 0.6. Because β controls the shape of the OCT intensity distribution, this means that it gets its most stretched shape in the IS/RPE and the least stretched in the ONL. According to our cellular hypothesis for the origin of stretched exponential, low β means that the IS/RPE expresses the most variation of reflectance across its local cellular-defined regions. Finally, we point out that although λ appears to be constant inside individual layers, it still shows substantial variation. Inspection of Fig. 6(a) reveals possible swings of up to 0.5 log units in λ within a layer. Future statistical studies may have to focus on these variations of λ to refine the ability to detect anomalies in diseased retinas.

D. Hope for Early Detection of Pathologies

The relative constancy of stretched exponential parameters, especially of k, within individual retinal layers of normal subjects gives hope that one may use OCT to detect pathologies automatically. One may even hope to detect pathologies with higher sensitivity and earlier than now possible. For such detection, certain pathologies should modulate these parameters by more than the variation observed in normal retinas. Hence, one may perhaps use statistical tools of detection. We are not the first to propose the use of statistical tools in the analysis of OCT data. Several groups have proposed statistical averaging-like tools to reduce speckle [55]–[59]. However, those methods are statistically suboptimal, since they focus on the first moments of the speckle intensity distribution. By fitting the probabilistic models developed in this paper to the retinal OCT data, we would be using all the information on their variability (i.e., on the speckle), not just the first moments. Using such information would bypass the necessity to smooth with arbitrary operators, leading to better algorithms.

For instance, the results of Fig. 7 suggest that such statistical tools may be useful to detect diabetic retinopathy. Modulations of λ and the normalization constant k of the stretched exponential distribution appear to provide enough signals to detect this disorder. For example, the SNR of the normalization constant k rises above 10 in Fig. 7(c). In contrast, diabetic retinopathy does not appear to affect the shape of the OCT intensity distribution, maintaining β relatively constant [see Fig. 7(b)].

How sensitive could a statistical method of detection be, given our findings of stretched exponential fits and of nonzero autocorrelation? To attempt answering this question, we performed pilot experiments with both real and simulated data. With real data, we tried to detect anomalies by exploiting the k curve within individual layers (Section II-D). We determined the baseline k and then considered anomalies to be ten times this value. With this simple method, samples in normal retina yielded less than 2% false-positive rate. In contrast, the hit rate (i.e., 100%—false-negative rate) for DME retinas was 93%, when detecting anomalies by eye. This percentage was obtained by viewing the pseudoimages before application of the algorithm and after extensive training. The 93% value applied to all anomalies in 34 images from each of eight DME patients for 272 pseudoimages. This high percentage and the low false-alarm rate were by no means definitive proof that the blots were markers of DME. However, these results suggested that DME did appear to produce spikes in the k curve. If confirmed, we will need further investigation to find out what these spikes are. They may mark zones of apoptosis, or may be due to exudates, hemorrhages, or microaneurysms in normal inner layers casting shadow on the outer layers [60]. Because one can distinguish apoptotic zones well with color fundus photography [61]–[63], we plan to perform a future study of these features combining this technique and OCT.

We also performed pilot studies with simulated “artificial retinas.” All pixels in the simulated control retinas had the same stretched exponential distribution intensities, with autocorrelation given as in the model of Section II-C1. Consequently, each simulated “sick retina” had, in addition, a randomly placed small region with lowered λ [see Fig. 7(a)]. This small region was our “artificial pathology.” “Normal” and “sick” retinas were generated at random, and we used a Bayesian method to detect the pathology. We then compared the level of Bayesian detection with that obtained by untrained human observers, using standard psychophysical techniques. Regardless of the β used in the experiment (0.5 ≤ β ≤ 1), Bayesian-detection thresholds (in terms of λ) were more than five times lower than human thresholds.

A possible explanation for the human disadvantage in detection threshold with artificial retinas is that statistics of natural images are different from those described in this paper [64]–[66]. The human visual system is adapted to natural statistics [64], [67], [68]. Consequently, human detection of stretched exponentials, which are common in natural images, may be impaired relative to well-designed statistical methods. Following our pilot results, an interesting question for future research is whether the main effect of the advance of the disease is a local, progressive lowering of λ and/or increasing of k of the stretched exponential fits. If so, because statistical methods are better than what one can achieve with the naked eye, then we may hope for improvement in detection even when the disease is in its early stages.

As a warning, we emphasize that this so-called human disadvantage is too premature for us to take it at face value. It refers only to psychophysical detection thresholds using “artificial retinas.” We must do much more work to compare our algorithms to human-expert performance using real OCT data. For instance, the methodology for human comparison would have to specify how observers were trained, how many of them were used, and how often images were observed. We plan to use such methodologies, but beforehand, we must develop optimal algorithms. The point of this paper is that before one develops optimal anomaly detection algorithms, one must develop good statistical models for how the normal data behave. Our stretched exponential model with its associated autocorrelation structure is exactly that. Because a nonoptimal algorithm based on this model outperformed human observers in a pilot, simplistic test is no proof that we are ready. However, this performance tantalizingly suggests that such statistical algorithms are worth pursuing.

ACKNOWLEDGMENT

The work of N. M. Grzywacz was supported in part by the National Eye Institute under Grant EY11170 and Grant EY016093 and in part by the National Science Foundation under Grant EEC-0310723. The work of D. Huang was supported by the National Eye Institute under Grant EY013516 and Grant P30 EY03040, and by a Grant from the Research to Prevent Blindness.

The authors are indebted to F. Spizzichino for pointing out the work of Langberg et al.. They also would like to thank S. Chatterjee, E. J. Lee, X. Cao, J. Lee, and J. Rapela for comments during performance of the work, and D. Steiner for administrative assistance. They performed part of this work during N. M. Grzywacz’s sabbatical stay at the Dipartimento di Informatica e Sistematica “A.Ruberti,” Università di Roma “La Sapienza.” N. M. Grzywacz wishes to thank the department for both its support and its hospitality.

Contributor Information

Norberto Mauricio Grzywacz, Departments of Biomedical and Electrical Engineering, Center for Vision Science and Technology, and the Neuroscience Graduate Program, University of Southern California, Los Angeles, CA 90089 USA (ude.csu.rsmb@gmn)

Joaquín de Juan, Departamento de Biotecnología, Universidad de Alicante, E-03080 Alicante, Spain (se.au@jdj)

Claudia Ferrone, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; ti.1amorinu.sid@enorref..

Daniela Giannini, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; ti.ilacsit@ininnaig.aleinad.

David Huang, Doheny Eye Institute and the Department of Ophthalmology, University of Southern California, Los Angeles, CA 90033 USA ; gro.ynehod@gnauhd.

Giorgio Koch, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; ti.orebil@oznerol_nas.

Valentina Russo, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; ti.1amorinu.sid@ossurv.

Ou Tan, Doheny Eye Institute and the Department of Ophthalmology, University of Southern California, Los Angeles, CA 90033 USA ; gro.ynehod@nato.

Carlo Bruni, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; ti.1amorinu.sid@inurb.olrac.

REFERENCES

1. Shapley R. Visual sensitivity and parallel retinocortical channels. Annu. Rev. Psychol. 1990;vol. 41:635–658. [PubMed]
2. Lee BB. Single units and sensation: A retrospect. Perception. 1999;vol. 28:1493–1508. [PubMed]
3. Grzywacz NM, Balboa RM. A Bayesian framework for sensory adaptation. Neural Comput. 2002 Mar.vol. 14:543–559. [PubMed]
4. Sterling P. Needle from a haystack. Optimal signaling by a nonlinear synapse. Neuron. 2002 May;vol. 34:670–672. [PubMed]
5. Fong DS, Aiello LP, Ferris FL, III, Klein R. Diabetic retinopathy. Diabetes Care. 2004 Oct.vol. 27:2540–2553. [PubMed]
6. Frank RN. Diabetic retinopathy. N. Engl. J. Med. 2004 Jan.vol. 350:48–58. [PubMed]
7. Jager RD, Mieler WF, Miller JW. Age-related macular degeneration. N. Engl. J. Med. 2008 Jun.vol. 358:2606–2617. [PubMed]
8. Sancho-Pelluz J, Arango-Gonzalez B, Kustermann S, Romero FJ, van Veen T, Zrenner E, Ekstrom P, Paquet-Durand F. Photoreceptor cell death mechanisms in inherited retinal degeneration. Mol. Neurobiol. 2008 Dec.vol. 38:253–269. [PubMed]
9. Resnikoff S, Pascolini D, Etya’ale DD, Kocur I, Pararajasegaram R, Pokharel GP, Mariotti SP. Global data on visual impairment in the year 2002. Bull. World Health Organ. 2004 Nov.vol. 82:844–851. [PubMed]
10. Sinclair SH. Diabetic retinopathy: The unmet needs for screening and a review of potential solutions. Expert Rev. Med. Devices. 2006 May;vol. 3:301–313. [PubMed]
11. Fong DS, Aiello L, Gardner TW, King GL, Blankenship G, Cavallerano JD, Ferris FL, III, Klein R. Retinopathy in diabetes. Diabetes Care. 2004 Jan.vol. 27:S84–S87. [PubMed]
12. Munier A, Gunning T, Kenny D, O’Keefe M. Causes of blindness in the adult population of the Republic of Ireland. Br. J. Ophthalmol. 1998 Jun.vol. 82:630–633. [PMC free article] [PubMed]
13. Harman-Boehm I, Sosna T, Lund-Andersen H, Porta M. The eyes in diabetes and diabetes through the eyes. Diabetes Res. Clin. Practice. 2007;vol. 78S:S51–S58.
14. Harris MI. Undiagnosed NIDDM: Clinical and public health issues. Diabetes Care. 1993 Apr.vol. 16:642–652. [PubMed]
15. Klein R. Prevention of visual loss from diabetic retinopathy. Surv. Ophthalmol. 2002 Dec.vol. 47:S246–S252. [PubMed]
16. Witkin SR, Klein R. Ophthalmologic care for persons with diabetes. J. Amer. Med. Assoc. 1984 May;vol. 251:2534–2537. [PubMed]
17. Costa RA, Skaf M, Melo LA, Jr, Calucci D, Cardillo JA, Castro JC, Huang D, Wojtkowski M. Retinal assessment using optical coherence tomography. Prog. Retinal Eye Res. 2006 May;vol. 25:325–353. [PubMed]
18. Sadda SR, Tan O, Walsh AC, Schuman JS, Varma R, Huang D. Automated detection of clinically significant macular edema by grid scanning optical coherence tomography. Ophthalmology. 2006 Jul.vol. 113:1187 e1–1187 e12. [PMC free article] [PubMed]
19. Virgili G, Menchini F, Dimastrogiovanni AF, Rapizzi E, Menchini U, Bandello F, Chiodini RG. Optical coherence tomography versus stereoscopic fundus photography or biomicroscopy for diagnosing diabetic macular edema: A systematic review. Invest. Ophthalmol. Vis. Sci. 2007 Nov.vol. 48:4963–4973. [PubMed]
20. Drexler W, Fujimoto JG. State-of-the-art retinal optical coherence tomography. Prog. Retinal Eye Res. 2008 Jan.vol. 27:45–88. [PubMed]
21. Dowling JE. Organization of vertebrate retinas. Invest. Ophthalmol. 1970 Sep.vol. 9:655–680. [PubMed]
22. Dacey DM. Primate retina: Cell types, circuits and color opponency. Prog. Retinal Eye Res. 1999 Nov.vol. 18:737–763. [PubMed]
23. Masland RH. Neuronal diversity in the retina. Curr. Opin. Neurobiol. 2001 Aug.vol. 11:431–436. [PubMed]
24. Goodman JW. Some fundamental properties of speckle. J. Opt. Soc. Amer. 1976;vol. 66:1145–1150.
25. Miller PF. The probability distribution of a wave at a very large depth within an extended random medium. J. Phys. A, Math. Gen. 1978;vol. 11:403–422.
26. Parry G. Measurement of atmospheric turbulence induced intensity fluctuations in a laser beam. Opt. Acta. 1981;vol. 28:715–728.
27. Schmitt JM, Xiang SH, Yung KM. Speckle in optical coherence tomography. J. Biomed. Opt. 1999;vol. 4:95–105. [PubMed]
28. Barber GW. Physiological chemistry of the eye. Arch. Ophthalmol. 1972 Jan.vol. 87:72–106. [PubMed]
29. Bishop PN. Structural macromolecules and supramolecular organisation of the vitreous gel. Prog. Retinal Eye Res. 2000 May;vol. 19:323–344. [PubMed]
30. Meek KM, Boote C. The organization of collagen in the corneal stroma. Exp. Eye Res. 2004 Mar.vol. 78:503–512. [PubMed]
31. Le Goff MM, Bishop PN. Adult vitreous structure and postnatal changes. Eye. 2008 Oct.vol. 22:1214–1222. [PubMed]
32. Barlow RE, Proschan F. Distributions with monotone failure rate. In: Serra A, Barlow RE, editors. Theory of Reliability (Società Italiana di Fisica XCIV Corso) Varenna, Italy: North Holland; 1986.
33. Langberg NA, Leon RV, Lynch J, Proschan F. The extreme points of the set of decreasing failure rate distributions. Z. Wahrschein-lichkeitstheorie Verw Gebiete. 1981;vol. 57:303–310.
34. Berberan-Santos MN, Bodunov EN, Valeur B. Mathematical functions for the analysis of the luminescence decays with underlying distributions 1. Kohlrausch decay function (stretched exponential) Chem. Phys. 2005;vol. 315:171–182.
35. Dunn OJ, Clark VA. Applied Statistics: Analysis of Variance and Regression. New York: Wiley; 1987.
36. Gloesmann M, Hermann B, Schubert C, Sattmann H, Ahnelt PK, Drexler W. Histologic correlation of pig retina radial stratification with ultrahigh-resolution optical coherence tomography. Invest. Ophthalmol. Vis. Sci. 2003 Apr.vol. 44:1696–1703. [PubMed]
37. Kim BY, Smith SD, Kaiser PK. Optical coherence tomographic patterns of diabetic macular edema. Amer. J. Ophthalmol. 2006 Sep.vol. 142:405–412. [PubMed]
38. Fischler MA, Bolles RC. Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. Commun. ACM. 1981;vol. 24:381–395.
39. Weibull W. A statistical distribution function of wide applicability. J. Appl. Mech. 1951;vol. 18:293–297.
40. Provis JM. The distribution and size of ganglion cells in the regina of the pigmented rabbit: A quantitative analysis. J. Comp. Neurol. 1979 May;vol. 185:121–137. [PubMed]
41. Stone J, Johnston E. The topography of primate retina: A study of the human, bushbaby, and new- and old-world monkeys. J. Comp. Neurol. 1981 Feb.vol. 196:205–223. [PubMed]
42. Rowe MH, Dreher B. Retinal W-cell projections to the medial interlaminar nucleus in the cat: Implications for ganglion cell classification. J. Comp. Neurol. 1982 Jan.vol. 204:117–133. [PubMed]
43. Fukuda Y, Hsiao CF, Watanabe M, Ito H. Morphological correlates of physiologically identified Y-, X-, and W-cells in cat retina. J. Neurophysiol. 1984 Dec.vol. 52:999–1013. [PubMed]
44. Tan O, Li G, Lu AT, Varma R, Huang D. Mapping of macular substructures with optical coherence tomography for glaucoma diagnosis. Ophthalmology. 2008 Jun.vol. 115:949–956. [PMC free article] [PubMed]
45. Williams G, Watts DC. Non-symmetrical dielectric relaxation behavior arising from a simple empirical decay function. Trans. Faraday Soc. 1970;vol. 66:80–85.
46. Laherrère J, Sornette D. Stretched exponential distributions in nature and economy: Fat tails with characteristic scales. Eur. Phys. J. B. 1998;vol. 2:525–539.
47. Dowling JE. Foveal receptors of the monkey retina: Fine structure. Science. 1965 Jan.vol. 147:57–59. [PubMed]
48. Braekevelt CR. Retinal photoreceptor fine structure in the domestic sheep. Acta Anat. Suppl. (Basel) 1983;vol. 116:265–275. [PubMed]
49. Yuodelis C, Hendrickson A. A qualitative and quantitative analysis of the human fovea during development. Vis. Res. 1986;vol. 26:847–855. [PubMed]
50. Dowling JE, Werblin FS. Synaptic organization of the vertebrate retina. Vision Res. 1971:1–15. [PubMed]
51. Wassle H, Yamashita M, Greferath U, Grunert U, Muller F. The rod bipolar cell of the mammalian retina. Vis. Neurosci. 1991 Jul-Aug;vol. 7:99–112. [PubMed]
52. Newman E, Reichenbach A. The Muller cell: A functional element of the retina. Trends Neurosci. 1996 Aug.vol. 19:307–312. [PubMed]
53. Sterling P, Cohen E, Freed MA, Smith RG. Microcircuitry of the on-beta ganglion cell in daylight, twilight, and starlight. Neurosci. Res. Suppl. 1987;vol. 6:S269–S285. [PubMed]
54. Wassle H, Boycott BB. Functional architecture of the mammalian retina. Physiol. Rev. 1991 Apr.vol. 71:447–480. [PubMed]
55. Ozcan A, Bilenca A, Desjardins AE, Bouma BE, Tearney GJ. Speckle reduction in optical coherence tomography images using digital filtering. J. Opt. Soc. Amer. A, Opt Image Sci. Vis. 2007;vol. 24:1901–1910. [PMC free article] [PubMed]
56. Abd-Elmoniem KZ, Youssef AB, Kadah YM. Real-time speckle reduction and coherence enhancement in ultrasound imaging via nonlinear anisotropic diffusion. IEEE Trans. Biomed. Eng. 2002 Sep.vol. 49(no. 9):997–1014. [PubMed]
57. Adler DC, Ko TH, Fujimoto JG. Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter. Opt. Lett. 2004;vol. 29:2878–2880. [PubMed]
58. Bashkansky M, Reintjes J. Statistics and reduction of speckle in optical coherence tomography. Opt. Lett. 2000;vol. 25:545–547. [PubMed]
59. Kim J, Miller DT, Kim E, Oh S, Oh J, Milner TE. Optical coherence tomography speckle reduction by a partially spatially coherent source. J. Biomed. Opt. 2005;vol. 10:064034. [PubMed]
60. Garcia JP, Jr, Garcia PM, Buxton DE, Panarelli A, Rosen RB. Imaging through opaque corneas using anterior segment optical coherence tomography. Ophthalmic Surg. Lasers Imag. 2007 Jul-Aug;vol. 38:314–318. [PubMed]
61. Niemeijer M, van Ginneken B, Russell SR, Suttorp-Schulten MS, Abramoff MD. Automated detection and differentiation of drusen, exudates, and cotton-wool spots in digital color fundus photographs for diabetic retinopathy diagnosis. Invest. Ophthalmol. Vis. Sci. 2007 May;vol. 48:2260–2267. [PMC free article] [PubMed]
62. Walter T, Massin P, Erginay A, Ordonez R, Jeulin C, Klein JC. Automatic detection of microaneurysms in color fundus images. Med. Image Anal. 2007 Dec.vol. 11:555–566. [PubMed]
63. Aptel F, Denis P, Rouberol F, Thivolet C. Screening of diabetic retinopathy: Effect of field number and mydriasis on sensitivity and specificity of digital fundus photography. Diabetes Metab. 2008 Jun.vol. 34:290–293. [PubMed]
64. Field DJ. Relations between the statistics of natural images and the response properties of cortical cells. J. Opt. Soc. Amer. A. 1987 Dec.vol. 4:2379–2394. [PubMed]
65. Ruderman DL, Bialek W. Statistics of natural images: Scaling in the woods. Phys. Rev. Lett. 1994 Aug.vol. 73:814–817. [PubMed]
66. Balboa RM, Grzywacz NM. Power spectra and distribution of contrasts of natural images from different habitats. Vis. Res. 2003 Nov.vol. 43:2527–2537. [PubMed]
67. Balboa RM, Grzywacz NM. The minimal local-asperity hypothesis of early retinal lateral inhibition. Neural Comput. 2000 Jul.vol. 12:1485–1517. [PubMed]
68. Simoncelli EP, Olshausen BA. Natural image statistics and neural representation. Annu. Rev. Neurosci. 2001;vol. 24:1193–1216. [PubMed]