|Home | About | Journals | Submit | Contact Us | Français|
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.
The retina is the first station of the brain’s visual system, and thus important for vision –. 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 , , age-related macular degeneration , and retinitis pigmentosa . In particular, diabetic retinopathy and macular degeneration are reaching epidemic levels . 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 . To emphasize the importance of early diagnosis, we consider the example of diabetic retinopathy. This disease is the most common microvascular complication of diabetes  and one of the leading causes of blindness in the western world . At the time of diagnosis, as many as 38% of type 2 diabetics manifest some degree of diabetic retinopathy . Diabetic retinopathy may precede the diagnosis of diabetes by as many as seven years and progress with the duration of diabetes . 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 . 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 . However, patients are not being diagnosed and treated in a timely fashion .
Optical coherence tomography (OCT) has recently become one of the primary methods for noninvasive mapping of the human retina –. Such maps can be especially useful for retinal analysis, as the retina has a strict layer organization –. 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 –. 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 –. In addition, some retinal layers contain large cell bodies of different types –, 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.
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.
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).
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
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.
Sufficiently small pixels of the pseudoimage should give rise to an exponential distribution of intensities –. 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
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
Hence, in particular, the mean, second moment, variance, and coefficient of variation are
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
In particular, for a stretched exponential distribution, the failure rate is
which is nonincreasing with respect to x [0,∞). Therefore, this distribution belongs to the class of the decreasing-failurerate (DFR) distributions . In turn, DFR distributions are intimately connected to the class of mixture-of-exponentials distributions . The following theoretical results hold for these distributions , .
A mixture of exponential distributions, i.e., a distribution characterized by the density
where µ is a probability measure on [0,∞), is DFR.
Moreover, on the basis of results proved by Langberg et al. , the following statement holds.
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
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:
Hence, using (4), one gets
Berberan–Santos et al.  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.
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  (p < 0.001).
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.,
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
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 , ρ, and ω, such that
Then, given two pairs of nonnegative numbers λ1 ≠ λ2 and β1 ≠ β2, we define the pairs (λx,βx) and (λy,βy) as follows:
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:
where δ denotes the Dirac delta function. One can show that µ is a probability measure provided that
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),µ(τ1,τ2) = µx(τ1)µy(τ2), 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 (x,y) as in (15), with marginals px(x) and py(y). From these functions, we define the joint density
As with µ, one can show that p(x,y) is a probability density with marginals px(x) and py(y), provided that
The parameter controls the correlation at the level of the measurements themselves. If = 0, then from (22) and the definition of (x,y), p(x,y) arises directly from µ through (15). The larger gets, the more the δ function dominates (22), making the occurrence of x = y more likely. Consequently, we can interpret 1 – 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
To calculate E(xy), we start from (22) to write
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
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
i.e., as expected, if the probability that the two measurements are equal goes to 1, so does the correlation. Fourth, if → 0, then R becomes proportional to ρ, and if ρ → 0, then R → . 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
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) = , because Γ(3)Γ(1) = 2Γ2(2) Thus, in particular
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 . We expect that, as the magnitude of 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 control correlation, [(28) or (30)], we propose that they fall with ||. For example, we can rewrite (30) as
This fall is consistent with the interpretation of the parameters given after (21) and (23). In particular, ρ should decrease with ||, as the measurement would be made from different cellular processes. Similarly, should decrease with ||, 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 to fall with ||. According to these physical interpretations for and ρ, we propose that the former falls faster with || 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 and ρ as the exponential decays
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 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 || > 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).
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.
To find these four or five borders, all algorithms broadly measured zones of statistically constant reflectance and then detected sharp gradients between these zones.
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 . 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).
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.
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 . 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).
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.
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.
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 , and their associated parameters υρ and υ. 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 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 υ are the lengths constant in our model, this expectation imposes υρ > υ. As Fig. 4(c) and (d) shows, fitting the data by optimizing υ and constraining υρ/υ = 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 –. 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 υρ/υ 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 yielded underestimation of correlation. For instance, for Fig. 4(c), when = 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 alone could yield better estimates than ρ alone, using only 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.
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.
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)].
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.
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 . However, blots are specially amenable to the techniques we developed here. (We plan to expand efforts  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.
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.
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 , , . 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 –. They cause the phase of individual rays of light to undergo a random walk –. 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.
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 , , . 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 – (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 = 0), then R = 0 (see also Section II-C2). (The reason for the necessity of the condition =0 is that ≠ 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 –. 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)].
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 , , . Other studies proposed a cellular identification of OCT layers , , , 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.
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 –. 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 . Because one can distinguish apoptotic zones well with color fundus photography –, 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 –. The human visual system is adapted to natural statistics , , . 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.
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.
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 (Email: ude.csu.rsmb@gmn)
Joaquín de Juan, Departamento de Biotecnología, Universidad de Alicante, E-03080 Alicante, Spain (Email: se.au@jdj)
Claudia Ferrone, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; Email: ti.1amorinu.sid@enorref..
Daniela Giannini, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; Email: email@example.com.
David Huang, Doheny Eye Institute and the Department of Ophthalmology, University of Southern California, Los Angeles, CA 90033 USA ; Email: gro.ynehod@gnauhd.
Giorgio Koch, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; Email: ti.orebil@oznerol_nas.
Valentina Russo, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; Email: ti.1amorinu.sid@ossurv.
Ou Tan, Doheny Eye Institute and the Department of Ophthalmology, University of Southern California, Los Angeles, CA 90033 USA ; Email: gro.ynehod@nato.
Carlo Bruni, Dipartimento di Informatica e Sistematica “A. Ruberti, ” Università di Roma “La Sapienza,” 00185 Rome, Italy ; Email: firstname.lastname@example.org.