Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Statistical analysis of corneal OCT speckle: a non-parametric approach

Open Access Open Access

Abstract

In biomedical optics, it is often of interest to statistically model the amplitude of the speckle using some distributional approximations with their parameters acting as biomarkers. In this paper, a paradigm shift is being advocated in which non-parametric approaches are used. Specifically, a range of distances, evaluated in different domains, between an empirical non-parametric distribution of the normalized speckle amplitude sample and the benchmark Rayleigh distribution, is considered. Using OCT images from phantoms, two ex-vivo experiments with porcine corneas and an in-vivo experiment with human corneas, an evidence is provided that the non-parametric approach, despite its simplicity, could lead to equivalent or better results than the parametric approaches with distributional approximations. Concluding, in practice, the non-parametric approach should be considered as the first choice to speckle modeling before a particular distributional approximation is utilized.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Speckle in optical coherence tomography (OCT) is usually treated as a source of intrinsic interference that needs to be reduced or eliminated [1,2]. Recently, a substantial interest in treating the OCT speckle as the source of information has been observed [37]. In such an approach, a particular local part of an OCT pseudo B-scan corresponding to the speckle amplitude signal, often termed a region of interest or a region of analysis, is treated as a random variable when the first-order statistics are considered. There is a convincing theoretical evidence, when employing the central limit theorem, that when the number of scattering elements in an observed sample is sufficiently large, such a random variable should be Rayleigh distributed, whereas when that number of scatterers is small, its statistics could be well approximated by the K distribution [4]. In fact, if the number of scatterers in a sample volume exceeds several elements, this already produce a distribution close to Rayleigh [8]. Despite that, a considerable effort has been made to model the OCT speckle amplitude with other types of distributions, often off-the-shelf and physically less justified. Examples of such, in an alphabetical order, include: Bessel K form [9], the physically justified special case of the Burr type XII distribution with one of its shape parameters, $c$, set to two [10], gamma [11], generalized extreme value [12], generalized gamma [5,13], log-normal [14], Nakagami [5], Rician [15], and Weibull [12,16] distributions. The main justification for using such models are their goodness-of-fit properties, which, in particular application, may achieve superior fitting performance to that of the two theoretically justified fundamental models [5,10,13], i.e., the Rayleigh and K distributions.

The purpose of modeling statistical properties of OCT speckle is to characterize its random amplitude with a set of distributional parameters that might be used for a particular practical task such as, for example, discriminating between different tissue samples. In that sense, in biomedical optics, the distributional parameters could then be viewed as biomarkers. However, when faced with the choice of an optimal (in some sense) distributional model, focus seems to be mainly given to the goodness-of-fit properties of a given model rather than to its physical correctness or the statistical power of model parameters that should subsequently act as biomarkers. Such an approach favors distributions with a larger number of parameters than one, even in cases when information criteria, such as the Akaike Information Criterion [5,10], are employed to account for over-parameterization. Further, in order to fit the amplitude of the OCT speckle with a particular distributional model, one needs to use a parameter estimation procedure, often based on the method of maximum likelihood, leading in case of distributions with more than one parameter to correlated parameter estimates, because orthogonality of the distributional parameters is rarely maintained. In fact, constructing distributions with more than two orthogonal parameters is not practical because the complexity of the differential equations involved in such a construction is high [17]. Hence, fitting of a three-parameter distribution, such as the Bessel K form [9] or the generalized gamma distribution [5,13], to the random OCT speckle amplitude can lead to a set of non-unique solutions, because different sets of parameter estimates could result in similarly well-fitted distributional representations. However, it should be emphasized that a better goodness-of-fit achieved with a particular distribution does not necessarily correspond to a better discriminating power of its parameter estimators [6].

In this paper, a paradigm shift is advocated to deviate from the parametric distributional approach to modeling the amplitude of the OCT speckle. Instead, after speckle amplitude normalization, the non-parametric approach based on the distance between the Rayleigh distribution with a fixed parameter $\sigma$ =$\sqrt {2}/2$ and an empirical distribution of the sample is considered. Such distances are being evaluated in the domain of probability density function, cumulative distribution function as well as that of the characteristic function. Also, the difference (distance) between the sample contrast ratio and the theoretical contrast ratio value for the Rayleigh distribution, is considered. To support our developments, OCT images from purposely designed phantoms, ex-vivo examination of porcine corneas and in-vivo examination of human corneas, are utilized. Further, the advocated non-parametric approach is contrasted against the parametric approach utilizing several distributional models.

2. Methods

Before introducing the incorporated image analysis methods involved in the proposed distributional model-free approach to speckle modeling, a summary of speckle theory and the associated statistical tools is first provided, for completeness.

2.1 Preliminaries

The fundamental theory of a speckle pattern formation was described by Goodman [18,19]. It is assumed that the electric field amplitude at some observation point of a speckle pattern results from contributions from different regions of the scattering media. The resultant complex amplitude is considered as a superposition of the elementary phasors with their amplitudes and phases statistically independent from each other and from amplitudes and phases of other phasors. The phases are also expected to be uniformly distributed $U(-\pi ,\pi )$. Having these assumptions, the real and imaginary parts of the resultant amplitude are considered to have zero mean, equal variances, and be uncorrelated. So if the number of elementary phasor contributions is large, the real and imaginary parts of the resultant phasor follow the Gaussian distribution, in compliance with the central limit theorem. It can be proven, that the length of the two-dimensional vector, which components are normally distributed with zero means and equal variances, $\sigma ^{2}$, obeys Rayleigh statistics. Therefore, the amplitude (length) $A$ of the resultant phasor in some point of the speckle field follows the Rayleigh distribution with the scale parameter $\sigma$, described by the probability density function (PDF) given by

$$p_A (A) = \frac{A}{\sigma ^{2}}e^{-\frac{A^{2}}{2\sigma ^{2}}}\:,\quad A\geq 0.$$

The aforementioned result is approximate for a finite large number of elementary phasors and becomes valid when the number of elementary phasors tends to infinity. The speckle field is then referred to as fully developed [20]. Additionally, for such a speckle field, the contrast ratio, defined as the ratio of the standard deviation of the amplitude to its mean value, approaches its Rayleigh-limited value of $\sqrt {4/\pi -1}\approx 0.5227$ [21].

The calculations of speckle statistics in OCT images are usually based on pixel values that represent the speckle field amplitude $A$ in a particular point (dependent on the image resolution) [22]. Since the speckle field amplitude follows the Rayleigh distribution, it is of interest to estimate the parameter of this distribution for the pixel values. In practice, it is usually considered to normalize the amplitude by dividing it by its root mean square (RMS) value [10,23], i.e., $A/\sqrt {\langle A^{2} \rangle }$. Such normalization, in many cases, leads to simplified parameter estimation procedures, in which one of the parameters (i.e., that of scale) achieves a particular constant value. It is worth noting that estimating the scale parameter for the normalized speckle amplitude may lead to incorrect interpretation of the estimated value because that estimator is biased and that bias depends on the other distributional parameter estimates (i.e., those of shape).

There are several methods available for estimating the scale parameter of the Rayleigh distribution. Because of its properties, the maximum likelihood estimator (MLE) is most frequently used and it takes the form

$$\hat{\sigma}_{\mathrm{MLE}}=\sqrt{\frac{1}{2n}\sum_{i=1}^{n}A_i^{2}}\:,$$
where $A_i$, $i=1,2,\dots , n$, with $n$ being the number of samples, are discrete samples of the speckle amplitude random variable $A$. If the amplitude of the speckle pattern is normalized, the above-mentioned estimator reduces to
$$\hat{\sigma}_{\mathrm{MLE}}=\sqrt{\frac{1}{2n}\cdot \frac{1}{\frac{1}{n}\sum_{i=1}^{n}A_i^{2}}\cdot \sum_{i=1}^{n}A_i^{2}}=\frac{\sqrt{2}}{2}\:,$$
leading to the Rayleigh distribution with a fixed scale parameter that is independent of the underlying distribution of the sample. Note that if the speckle was fully developed, its normalized amplitude would follow the Rayleigh distribution with a fixed scale parameter, but the converse is not true. One may consider another estimator for $\sigma$ than MLE such as that proposed by Ardianti who used a Bayes method of the form [24]
$$\hat{\sigma}_{\mathrm{Bayes}} = \frac{\sqrt{2}\Gamma(n+2)}{2\Gamma\big(n+\frac{5}{2}\big)}\sqrt{\sum_{i=1}^{n}A_i^{2}}\:.$$

Again, for the normalized speckle amplitude one obtains

$$\hat{\sigma}_{\mathrm{Bayes}} = \frac{\sqrt{2}\Gamma(n+2)}{2\Gamma\big(n+\frac{5}{2}\big)}\sqrt{\frac{1}{\frac{1}{n}\sum_{i=1}^{n}A_i^{2}}\cdot \sum_{i=1}^{n}A_i^{2}}=\frac{\sqrt{2}}{2}\cdot \frac{\Gamma(n+2)}{\Gamma\big(n+\frac{5}{2}\big)}\cdot \sqrt{n} \:,$$
which asymptotically, for $n\to \infty$, approaches $\sqrt {2}/2$. This indicates that for any random sample that is normalized, the estimate of $\sigma$ has a fixed value, regardless of the estimation method. Therefore, the Rayleigh distribution with a scale parameter equal to $\sqrt {2}/2$ can be used as a benchmark.

The distribution described above is theoretically justified for the speckle amplitude provided that there is a sufficient number of scatterers in a sample volume, defined by the OCT beam cross section and coherence length [8]. However, sometimes the number of scatterers may be insufficient to obtain a fully developed speckle field. The speckle statistics are then altered, as the central limit theorem is no longer applicable. The theoretical derivations on that problem, developed by Jakeman and Pusey [25], included the assumption that the number of contributions to the scattered field is fluctuating and could be modeled by a negative binomial distribution. That leads to the K distribution as a model of the speckle field amplitude with the form [26]

$$p_A(A)=\frac{4}{\Gamma (\alpha)}\sqrt{\frac{\alpha}{\langle A^2\rangle}}\cdot \Bigg( \frac{\alpha A^2}{\langle A^2\rangle} \Bigg)^\frac{\alpha}{2}K_{\alpha-1}\Bigg(2\sqrt{\frac{\alpha A^2}{\langle A^2\rangle}} \Bigg) \:,$$
where $\alpha$ is the shape parameter of the K distribution. This model is widely used in works where a small number of scatterers is considered. Also, the shape parameter of the K distribution is linked to the average scatterer density in the tested object [2729]. Owing to the fact that the Rayleigh distribution is theoretically proven to be suitable for modeling the fully developed speckle field amplitude, in section 2.4 we advocate a non-parametric approach for speckle statistics analysis, which does not involve any distributional parameter estimation. However, first, the experimental data, used to illustrate the proposed approach and its potential applications, is described.

2.2 Experimental data

To validate the non-parametric approach, data from three studies performed on resin phantoms, ex-vivo on porcine eyeballs, and in-vivo on human eyes, are used. For each of them, the OCT B-scans were acquired using a spectral OCT (SOCT Copernicus REVO, Optopol, Zawiercie, Poland), with the center wavelength of 830 nm, the half bandwidth of 50 nm, the axial resolution of 5 µm, and the transversal resolution of 15 µm. The scanning speed of the device is 80 000 A-scans per second. The measurements were acquired at a constant aperture within the bands of the instrument’s depth of focus using the own guiding system of the instrument.

2.2.1 Phantoms

The first study, performed on phantoms, was a pilot study to assess the impact of scatterer density on the speckle statistics and the results of the non-parametric approach. For this, a set of purposely designed phantoms were fabricated. They were made of the epoxy resin L-285 (Havel composites, Cieszyn, Poland), which in its liquid state was thoroughly mixed with the blue dye powder particles of approximate size of a few to tens micrometers. When the dye powder was visually uniformly distributed in the resin, the mixture was carefully poured on a microscope slide, to form a shape of a convex disc with the thickness of about 1 mm and the diameter of about 10 mm. After drying, the phantoms became solid transparent discs with visible blue particles. This procedure was repeated 9 times, increasing the amount of the powder dye, to obtain phantoms with different scattering particle concentrations (C1, C2,…, C9). The procedure was such that at C1 about 10 mm$^{3}$ of dye was added to 1 cm$^{3}$ of epoxy raisin. At C2, about 20 mm$^{3}$ was additionally added to it, then at C3 about 40 mm$^{3}$ was added, and so on, every time adding twice as much dye as in the previous step. The discs’ dimensions and the concentration of the dye were not precisely controlled but the increasing trend of concentration between subsequent phantoms was well maintained. For each phantom, a single OCT B-scan of size $3077\times 708$ pixels, encompassing the central part of the phantom, was registered.

2.2.2 Porcine corneas: ex-vivo study

Data from a recent study [7] consisting of two ex-vivo experiments on porcine eyeballs, in which the parametric distributional-based approach to speckle modeling was used, were taken for validating the non-parametric approach. Experiment 1 was performed to evaluate the influence of intraocular pressure (IOP) elevation on the corneal OCT speckle statistics. The IOP in the anterior chamber of each eyeball was increased from 10 mmHg to 40 mmHg with the step of 5 mmHg and the OCT scans of size $1536\times 736$ pixels were registered at each IOP level. Experiment 2 was prepared according to the same procedure, but with IOP maintained at the constant level of 15 mmHg. The OCT scans were collected in 7 time points, adequate to the first experiment, to examine the impact of the experiment’s duration on speckle statistics. Thirty three eyeballs were included in the study: 23 eyeballs for Experiment 1 and 10 eyeballs for Experiment 2. The entire procedure of the experiments’ preparation and the characteristics of the equipment are described in detail in that recent work [7].

2.2.3 Human corneas: in-vivo study

The final set of data, used here for the validation of the non-parametric approach to speckle analysis, were measurements of the central cornea in a group of 56 healthy Caucasian subjects with the mean ($\pm$ standard deviation) age of $42.4 \pm 18.3$ years (range from 21 to 87 years). For each subject, a single OCT B-scan of size $1536\times 736$ was acquired for randomly chosen eye and a non-corrected IOP was measured using the noncontact tonometer Corvis Scheimpflug Technology (Corvis ST, OCULUS, Wetzlar, Germany). All registered values of IOP were within the normal physiological limits (i.e., less than or equal to 20 mmHg). This part of the study was performed in accordance with the tenets of the Declaration of Helsinki.

2.3 Image analysis

All calculations were performed in MATLAB (MathWorks, Inc. Natick, MA, USA). At first, the log-transformation, applied automatically to every OCT B-scan in the device software, was inversed. Pixel values were selected from the specified region of interest (ROI), dependent on the type of imaged object (Fig. 1). For phantoms, the ROI was placed on the left of the central reflection, 20 pixels below the top border of the phantom, and was of size $600\times 220$ pixels, encompassing an area of about 2 mm horizontally by 0.6 mm vertically. Note that for imaging phantoms, an instrument protocol with an external adapting lens was used, allowing the so-called wide scan. For porcine and human corneas the ROI was placed centrally in relation to the apex, 10 pixels below the Bowman’s layer and 10 pixels above the endothelium, and had the width of 600 pixels, which corresponds to the central 2 mm of the cornea. For pixel values from each ROI, treated as a random variable $X$, empirical cumulative distribution function (eCDF), kernel density estimator (KDE), empirical characteristic function (eCF), and the sample contrast ratio (CR), were calculated in accordance with the formulas:

  • • eCDF [30]
    $$\mathrm{eCDF}=\frac{1}{n}\sum_{i=1}^{n}\textbf{1}\{X_i\leq t\} \:,$$
    where $\textbf {1}\{X_i\leq t\}$ is the indicator of the event that the value of a random variable $X_i$ is less than or equal to $t$, whereas $n$ is the number of samples in ROI.
  • • KDE [31]
    $$\mathrm{KDE}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{h}K\Big(\frac{x-X_i}{h}\Big)\:,$$
    where $h$ is the bandwidth of the estimator and $K$ is the non-negative kernel function. The Gaussian type of kernel function was set and the bandwidth value was $h=(0.75n)^{(-1/5)}\hat {\sigma }_X$ with $\hat {\sigma }_X$ being the sample standard deviation. In KDE calculations it was necessary to correct a boundary effect, related to the left side of the interval of pixel values, which are nonnegative. For this purpose, KDE values for abscissa less than 0.05 were not taken into account in further calculations [32].
  • • eCF [33]
    $$\mathrm{eCF}=\frac{1}{n}\sum_{i=1}^{n}e^{jtX_i}, \quad j=\sqrt{-1}\:.$$
  • • CR [19]
    $$\mathrm{CR} = \frac{\hat{\sigma}_X}{\overline{X}}\:,$$
    where $\overline {X}$ denotes the sample mean.

 figure: Fig. 1.

Fig. 1. Illustrative OCT scans from each considered study with ROI marked with red lines. Cyan lines indicate phantom borders (top image) or epithelium, Bowman’s layer, and endothelium (the middle and bottom image).

Download Full Size | PDF

2.4 Non-parametric approach

In this work, a non-parametric approach for OCT corneal speckle analysis is advocated, in which the Rayleigh distribution with a scale parameter $\sigma =\sqrt {2}/2$ is used as a benchmark distribution (BD). The empirical distributions of pixel values from the specified ROI in OCT images are compared with this benchmark using the following four different statistical distances:

  • • Kolmogorov–Smirnov distance between eCDF and the cumulative distribution function (CDF) of the benchmark distribution [34]
    $$D_{\mathrm{KS}}=\sup_x|\mathrm{eCDF}(x)-\mathrm{CDF_{BD}}(x)|$$
    where $\sup _x$ is the supremum of the set of distances across all $x$ values.
  • • mean square error (MSE) distance, between KDE and the probability density function (PDF) of the benchmark distribution
    $$D_{\mathrm{MSE}} = \frac{1}{n}\sum_{i=1}^{n}(\mathrm{KDE}(x_i)-\mathrm{PDF_{BD}}(x_i))^{2}$$
  • • maximum mean discrepancy (MMD) distance, between eCF and the characteristic function (CF) of the benchmark distribution [35]
    $$D_{\mathrm{MMD}}= \left\lVert{\frac{1}{n}\sum_{i=1}^{n}\mathrm{eCF}(x_i)-\frac{1}{n}\sum_{i=1}^{n}\mathrm{CF_{BD}}(x_i)}\right\rVert$$
  • • CR distance, between CR calculated from pixel values and the theoretical value of 0.5227 for the Rayleigh distribution
    $$D_{\mathrm{CR}}=\mathrm{CR}-0.5227$$

2.5 Statistical analysis

In the ex-vivo study on porcine corneas, changes in statistical distance values were evaluated for two experiments. One way repeated measures analysis of variance (rmANOVA) was used to assess whether the values of considered statistical distances vary with IOP (Experiment 1) and with time (Experiment 2). Post-hoc analysis was involved to assess the statistical significance of changes in distance values between adjacent IOP levels or time points. To obtain the correlation coefficient between the statistical distance and IOP or time, linear mixed effect model fitting was applied, because for the same eyeball distance values for consecutive IOP levels or time points are correlated. In in-vivo study on human corneas, Pearson correlation coefficients were calculated between distance values and IOP. Differences between correlation coefficients were tested using the Fisher test.

3. Results

3.1 Phantoms

Measurements of the fabricated phantoms were used to assess whether the scatterer density has an impact on the statistical distance values.

To justify the ROI analyses it was examined first whether the sample mean and sample standard deviation substantially change from place to place. For that, a sliding window (ROI) was moved towards the centre of each scan in steps of 10 pixels providing 30 additional ROIs. Each shifted ROI mean and standard deviation was compared to those of the originally selected ROI. The further the window was shifted, the larger deviations from the original ROI were noticed, indicating that the assumption of stationarity and ergodicity of ROI is not strictly maintained. However, examining a kind of an effect size, the maximum change in the mean did not exceed 5.5% for the lowest concentration (C1) and 2.2% for the highest concentration (C9) whereas the maximum change in the standard deviation did not exceed 5% for the highest concentration (C9) and 2.5% for lowest concentration (C1). Consequently, the assumption of stationarity, at least in terms of stable mean and stable standard deviation, could be made.

Figure 2 shows that the decrease in scatterer concentration causes deterioration of the Rayleigh model fitting to experimental data, meaning that the empirical distribution becomes more divergent from the benchmark distribution. Hence, the Rayleigh distribution with a scale parameter $\sqrt {2}/2$ becomes more valid when the number of scatterers is larger, which is consistent with theoretical derivations on speckle field amplitude, described in section 2.1. Contrast ratio distance, as well as other distances, have similar, less than 0.2, values for higher dye concentrations (C9 – C6) and starts to increase from concentration C5 achieving highest values for concentration C1. It is worth noting that by reducing the size of ROI to a small local neighborhood [14], for a fully developed speckle one can achieve distance values approaching zero. Figure 3 shows the estimated distances as functions of ROI size for the phantom of concentration C9.

 figure: Fig. 2.

Fig. 2. Distance values for phantoms with decreasing concentration of blue dye particles. Illustrative images inserted at the bottom of the plot present the ROIs of the OCT images of the phantoms for the concentrations C8, C5, and C2.

Download Full Size | PDF

 figure: Fig. 3.

Fig. 3. Estimated distances as functions of ROI area for the phantom of concentration C9. S denotes the original ROI area of $600\times 220$ pixels.

Download Full Size | PDF

To facilitate the comparison between the proposed non-parametric approaches and parametric distributional approaches, a set of eight models, namely the two-parameter Burr distribution (with one of its shape parameters, $c$, set to 2, as it was defined in [10]), the three-parameter Burr distribution, gamma, generalized gamma, K, Nakagami, Rayleigh, and Weibull distributions, were used. For each considered phantom ROI sample, the respective distributional parameters, here denoted as scale and shape parameters, were estimated using the MLE methods, except for the case of the two-parameter Burr distribution for which a method based on the sample mean and sample median was utilized. The results of this analysis are shown in Fig. 4, where similar trends with respect to scatterer concentration to those shown in Fig. 2 are evident, particularly for the shape parameters of the two-parameter distributions.

 figure: Fig. 4.

Fig. 4. Values of parameters of the considered distributions for decreasing concentrations of dye particles in the phantom study.

Download Full Size | PDF

3.2 Porcine corneas: ex-vivo study

Table 1 (columns 3 and 4) shows the results of repeated measures ANOVA for both non-parametric and parametric approaches for Experiment 1 and Experiment 2 whereas Fig. 5 presents the boxplots of the non-parametric distance values for two ex-vivo experiments on porcine eyes. There are statistically significant differences between non-parametric distance values for different IOP levels (p1) and for different time points (p2). Post-hoc analysis was performed to assess statistical significance of distance values differences between adjacent IOP levels and consecutive time points. It revealed that in Experiment 2 there are statistically significant differences in three of considered distances values between $t_4$ and $t_5$ time points ($p = 0.024$, $p = 0.049$ and $p = 0.020$ for $D_{\mathrm {MSE}}$, $D_{\mathrm {KS}}$ and $D_{\mathrm {MMD}}$, respectively) and in $D_{\mathrm {KS}}$ values between $t_3$ and $t_4$ ($p = 0.043$). Statistically significant differences are also observed in $D_{\mathrm {CR}}$ values for pairs of adjacent time points from $t_1$ to $t_4$ ($p < 0.001$, $p = 0.030$, $p = 0.004$). In Experiment 1, there are statistically significant differences in values of all distances between adjacent IOP levels from 15 mmHg to 35 mmHg (consecutive p-values were: for $D_{\mathrm {MSE}}$ 0.049, $<0.001$, 0.001, 0.008; for $D_{\mathrm {KS}}$ 0.025, $<0.001$, 0.001, 0.002; for $D_{\mathrm {MMD}}$ 0.023, $<0.001$, 0.002, 0.003; for $D_{\mathrm {CR}}$ $<0.001$, $<0.001$, 0.002, 0.002). The correlation coefficients between statistical distances and IOP, calculated using linear mixed-effect model fitting, are 0.841, 0.864, 0.844 and 0.880 for $D_{\mathrm {MSE}}$, $D_{\mathrm {KS}}$ $D_{\mathrm {MMD}}$ and $D_{\mathrm {CR}}$, respectively.

 figure: Fig. 5.

Fig. 5. Boxplots of $D_{\mathrm {MSE}}$, $D_{\mathrm {KS}}$, $D_{\mathrm {MMD}}$, and $D_{\mathrm {CR}}$ values for Experiment 1 (examining the impact of IOP) and Experiment 2 (examining the impact of its duration) on porcine corneas (ex-vivo study). Asterisks indicate statistical significance of differences between the values of statistical distances for adjacent IOP levels (blue ones) or time points (red ones), assessed using paired t-test (* $p<0.05$; ** $p<0.01$; *** $p<0.001$).

Download Full Size | PDF

Tables Icon

Table 1. The results of repeated measures ANOVA for two ex-vivo porcine eye experiments and the results of correlation analysis for the in-vivo human eye study for both non-parametric and parametric approaches.

Since, as noted in the Introduction, a better goodness-of-fit achieved with a particular parametric distribution does not necessarily correspond to a better discriminating power of its parameter estimators, the considered set of eight distributions was used to fit the data from the first ex-vivo porcine eye study (Experiment 1). Further, for each measurement, the goodness-of-fit (GoF), here corresponding to the mean square error between the given estimated PDF and the KDE, was calculated. The results of the GoF are shown in Fig. 6. It is evident that the two three-parameter distributions, namely, the generalized gamma and the Burr, achieve better GoF results than those of one- or two-parameter distributions. Interestingly, the two-parameter Burr distribution achieves the best GoF among other one- and two-parameter distributions.

 figure: Fig. 6.

Fig. 6. Means and standard errors of the goodness-of-fit (GoF) for a set of different distributions for the OCT images of porcine eyeballs in the ex-vivo study in Experiment 1.

Download Full Size | PDF

3.3 Human corneas: in-vivo study

The in-vivo study on human corneas showed that there are statistically significant correlations between IOP and values of statistical distances (see columns 5 and 6 of Table 1, where the values of Pearson correlation coefficient as well as corresponding p-values for testing the statistical significance of observed correlations are shown). The correlations were weak but showed that the overall trend of statistical distance values increasing with IOP, found in the porcine corneas in the ex-vivo study, is preserved. It is worth noting that the maximum correlation coefficient value of 0.401, achieved for $D_{\mathrm {MSE}}$, is not statistically significantly different to the minimum correlation coefficient value of 0.364, achieved for $D_{\mathrm {CR}}$ (Fisher test, $p = 0.396$). To compare the results of parametric and non-parametric approaches, correlations between the parameters of considered distributions and IOP were also calculated. For parametric distributions, the maximum correlation between distributional parameters and IOP is achieved for the two-parameter Burr distribution. Note that the maximum achieved correlation value of $-0.492$ for the shape parameter of the two-parameter Burr distribution is not statistically significantly different to the maximum coefficient value of 0.401, achieved for $D_{\mathrm {MSE}}$ (Fisher test, $p = 0.055$) but it is statistically significantly different to the correlation coefficient values achieved for the remaining statistical distances (Fisher test, $p = 0.047$, $p = 0.034$ and $p = 0.022$ for $D_{\mathrm {KS}}$, $D_{\mathrm {MMD}}$ and $D_{\mathrm {CR}}$, respectively).

4. Discussion

In this paper, several non-parametric statistical distances based on PDF, CDF, and CF, were employed to quantify the similarity of the empirical distribution of the normalized speckle amplitude to the benchmark Rayleigh distribution with the scale parameter $\sigma = \sqrt {2}/2$. Additionally, a distance between the sample contrast ratio and that of the Rayleigh distribution was considered. The Rayleigh distribution is an appropriate physically justified model for speckle field amplitude provided that the number of scatterers is relatively large. Hence, it is concluded that increasing the number of scatterers in a tested object should tend the empirical distribution of the speckle amplitude to the Rayleigh distribution, minimizing the statistical distance between those two distributions. Such an approach has been used by Matveev et al. [36] and Demidov et al. [37], who used a Rayleigh distribution goodness-of-fitting to isolate low scattering biological structures and to distinguish them from noise. Here, this study additionally contributes to those developments, showing that statistical distance measures, independent of the domain that are calculated in, can be successfully utilized to distinguish between different samples as well as different sample conditions.

Regarding the contrast ratio, which value for a fully developed speckle field, that can be modeled by the Rayleigh distribution, is approximately equal to 0.52, theoretical and experimental results, described by Hillman et al. [21], show the negligible variation in the contrast ratio for a large number of scatterers and strong variation of its values for lower number of scatterers. The use of the contrast ratio for scatterer density characterization, especially for the speckle fields that are not fully developed, was proposed earlier [14,38]. The results of the phantom experiment considered in this study confirm those developments. The rationale behind using a statistical distance rather than a particular distributional model itself is evident when considering the results for the fabricated phantoms. It is clear that for low concentration levels (C1, C2,…, C5) the speckle pattern is not fully developed (the CR distance does not approach zero) whereas for higher concentration levels (C6, C7,…, C9) that speckle approaches that of a fully developed pattern.

In the ex-vivo study on porcine eyes the increase of statistical distance values was observed during the IOP elevation. This suggests that the increase in IOP causes the decrease in scatterer density and that is in agreement with the study of Wu et al. [39] who showed similar evolution of collagen microstructure, in particular the immediate loss of interlamellar gaps, with increasing IOP using nonlinear optical microscopy. Hence, the evaluation of statistical distances for the speckle pattern can provide some insight into the microstructure of the imaged tissue sample. Furthermore, the in-vivo study on human corneas confirms the positive trend of statistical distance values with increasing IOP that has been shown earlier in the porcine ex-vivo study.

Further, the non-parametric approaches have been contrasted against parametric approaches using a set of eight distributions. This comparison clearly indicates that the performance of the advocated here non-parametric statistical distances is better than that of the majority of parametric approaches, with the only exception of the two-parameter Burr distribution, for which their performance, at least for the data in the considered studies, is equivalent. It is also important to note that the non-parametric approaches possess generally higher computational efficiency than those of the parametric techniques, for which numerical MLE methods need to be employed.

Summarizing, when considering the first order statistics of the speckle amplitude, the distance between the sample contrast ratio and that of the Rayleigh distribution is the simplest, yet equally effective measure of speckle departure from a fully developed field than other statistical distances evaluated in different domains. This study emphasizes that such a non-parametric approach should be considered before any other, often more complicated, parametric approach is utilized.

Funding

Narodowe Centrum Nauki (2018/29/B/ST7/02451); Narodowa Agencja Wymiany Akademickiej (PPI/APM/2019/1/00085/DEC/1); InterDok — Interdisciplinary Doctoral Studies Projects at Wroclaw University of Science and Technology.

Acknowledgments

Authors wish to express their gratitude and thank Dr Monika Danielewska and Dr Małgorzata Kostyszak for their joint contributions to the three considered studies included in this work as well as for the fruitful discussions. Authors thank the anonymous reviewers for their excellent constructive criticism as well as supportive feedback that allowed better shaping this work.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef]  

2. Y. Zhao, K. K. Chu, W. J. Eldridge, E. T. Jelly, M. Crose, and A. Wax, “Real-time speckle reduction in optical coherence tomography using the dual window method,” Biomed. Opt. Express 9(2), 616–622 (2018). [CrossRef]  

3. M. Y. Kirillin, G. Farhat, E. A. Sergeeva, M. C. Kolios, and A. Vitkin, “Speckle statistics in OCT images: Monte Carlo simulations and experimental studies,” Opt. Lett. 39(12), 3472–3475 (2014). [CrossRef]  

4. M. Almasian, T. G. Van Leeuwen, and D. J. Faber, “OCT amplitude and speckle statistics of discrete random media,” Sci. Rep. 7(1), 14873 (2017). [CrossRef]  

5. D. A. Jesus and D. R. Iskander, “Assessment of corneal properties based on statistical modeling of OCT speckle,” Biomed. Opt. Express 8(1), 162–176 (2017). [CrossRef]  

6. D. R. Iskander, M. A. Kostyszak, D. A. Jesus, M. Majewska, M. E. Danielewska, and P. Krzyżanowska-Berkowska, “Assessing corneal speckle in optical coherence tomography: a new look at glaucomatous eyes,” Optom. Vis. Sci. 97(2), 62–67 (2020). [CrossRef]  

7. M. Niemczyk, M. E. Danielewska, M. A. Kostyszak, D. Lewandowski, and D. R. Iskander, “The effect of intraocular pressure elevation and related ocular biometry changes on corneal OCT speckle distribution in porcine eyes,” PLoS One 16(3), e0249213 (2021). [CrossRef]  

8. V. Y. Zaitsev, L. Matveev, A. Matveyev, G. Gelikonov, and V. Gelikonov, “A model for simulating speckle-pattern evolution based on close to reality procedures used in spectral-domain OCT,” Laser Phys. Lett. 11(10), 105601 (2014). [CrossRef]  

9. M. Samieinasab, Z. Amini, and H. Rabbani, “Multivariate statistical modeling of retinal optical coherence tomography,” IEEE Trans. on Med. Imaging 39(11), 3475–3487 (2020). [CrossRef]  

10. G. R. Ge, J. P. Rolland, and K. J. Parker, “Speckle statistics of biological tissues in optical coherence tomography,” Biomed. Opt. Express 12(7), 4179–4191 (2021). [CrossRef]  

11. A. A. Lindenmaier, L. Conroy, G. Farhat, R. S. DaCosta, C. Flueraru, and I. A. Vitkin, “Texture analysis of optical coherence tomography speckle for characterizing biological tissues in vivo,” Opt. Lett. 38(8), 1280–1282 (2013). [CrossRef]  

12. R. Shetty, M. Francis, R. Shroff, N. Pahuja, P. Khamar, M. Girrish, R. M. Nuijts, and A. S. Roy, “Corneal biomechanical changes and tissue remodeling after SMILE and LASIK,” Invest. Ophthalmol. Visual Sci. 58(13), 5703–5712 (2017). [CrossRef]  

13. M. E. Danielewska, A. Antonczyk, D. A. Jesus, M. M. Rogala, A. Błonska, M. Cwirko, Z. Kiełbowicz, and D. R. Iskander, “Corneal optical coherence tomography speckle in crosslinked and untreated rabbit eyes in response to elevated intraocular pressure,” Trans. Vis. Sci. Tech. 10(5), 2–12 (2021). [CrossRef]  

14. D. D. Duncan, S. J. Kirkpatrick, and R. K. Wang, “Statistics of local speckle contrast,” J. Opt. Soc. Am. A 25(1), 9–15 (2008). [CrossRef]  

15. K. H. Cheng, A. Mariampillai, K. K. Lee, B. Vuong, T. W. Luk, J. Ramjist, A. Curtis, H. Jakubovic, P. Kertes, M. Letarte, M. E. Faughnan, and V. X. Yang, “Histogram flow mapping with optical coherence tomography for in vivo skin angiography of hereditary hemorrhagic telangiectasia,” J. Biomed. Opt. 19(8), 086015 (2014). [CrossRef]  

16. V. Kajic, “Automated retinal layer segmentation and pre-apoptotic monitoring for three-dimensional optical coherence tomography,” Doctoral dissertation, Cardiff University (2011).

17. V. S. Huzurbazar, “Probability distributions and orthogonal parameters,” Math. Proc. Cambridge Philos. Soc. 46(2), 281–284 (1950). [CrossRef]  

18. J. W. Goodman, “Statistical properties of laser speckle patterns,” in Laser Speckle and Related Phenomena, J. C. Dainty, ed. (Springer, 1975), pp. 9–75.

19. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (SPIE, 2007).

20. N. George, C. R. Christensen, J. S. Bennett, and B. D. Guenther, “Speckle noise in displays,” J. Opt. Soc. Am. 66(11), 1282–1290 (1976). [CrossRef]  

21. T. R. Hillman, S. G. Adie, V. Seemann, J. J. Armstrong, S. L. Jacques, and D. D. Sampson, “Correlation of static speckle with sample properties in optical coherence tomography,” Opt. Lett. 31(2), 190–192 (2006). [CrossRef]  

22. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography — principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]  

23. T. K. Stanton, W.-J. Lee, and K. Baik, “Echo statistics associated with discrete scatterers: A tutorial on physics-based methods,” J. Acoust. Soc. Am. 144(6), 3124–3171 (2018). [CrossRef]  

24. F. Ardianti, “Estimating parameter of Rayleigh distribution by using maximum likelihood method and Bayes method,” in IOP Conference Series: Materials Science and Engineering, vol. 300 (2018), p. 012036.

25. E. Jakeman and P. N. Pusey, “Significance of K distributions in scattering experiments,” Phys. Rev. Lett. 40(9), 546–550 (1978). [CrossRef]  

26. E. Jakeman, “Speckle statistics with a small number of scatterers,” Opt. Eng. 23(4), 453–461 (1984). [CrossRef]  

27. I. A. Popov, N. V. Sidorovsky, and L. M. Veselov, “Experimental study of intensity probability density function in the speckle pattern formed by a small number of scatterers,” Opt. Commun. 97(5-6), 304–306 (1993). [CrossRef]  

28. M. Sugita, A. Weatherbee, K. Bizheva, I. Popov, and A. Vitkin, “Analysis of scattering statistics and governing distribution functions in optical coherence tomography,” Biomed. Opt. Express 7(7), 2551–2564 (2016). [CrossRef]  

29. A. Weatherbee, M. Sugita, K. Bizheva, I. Popov, and A. Vitkin, “Probability density function formalism for optical coherence tomography signal analysis: a controlled phantom study,” Opt. Lett. 41(12), 2727 (2016). [CrossRef]  

30. A. W. van der Vaart, “Asymptotic statistics,” in Asymptotic Statistics (Cambridge University, 1998), chap. Empirical, pp. 265–290.

31. B. W. Silverman, Density Estimation for Statistics and Data Analysis (Routledge, 2018).

32. M. C. Jones, “Simple boundary correction for kernel density estimation,” Stat. Comput. 3(3), 135–146 (1993). [CrossRef]  

33. H. Cramer, Mathematical Methods of Statistics (Asia Publishing House, 1946).

34. F. J. Massey, “The Kolmogorov–Smirnov test for goodness of fit,” J. Am. Stat. Assoc. 46(253), 68–78 (1951). [CrossRef]  

35. A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola, “A kernel two-sample test,” J. Mach. Learn. Res. 13, 723–773 (2012).

36. L. A. Matveev, A. L. Matveyev, V. Demidov, A. A. Sovetsky, G. V. Gelikonov, V. Y. Zaitsev, and I. A. Vitkin, “Assessment of optical coherence tomography speckle patterns in low-scatterer-concentration regions: simulations for lymphatic vessels mapping,” in European Conference on Biomedical Optics (Optical Society of America, 2019), p. 11075_2.

37. V. Demidov, L. A. Matveev, O. Demidova, A. L. Matveyev, V. Y. Zaitsev, C. Flueraru, and I. A. Vitkin, “Analysis of low-scattering regions in optical coherence tomography: applications to neurography and lymphangiography,” Biomed. Opt. Express 10(8), 4207–4219 (2019). [CrossRef]  

38. G. Lamouche, B. F. Kennedy, K. M. Kennedy, C.-E. Bisaillon, A. Curatolo, G. Campbell, V. Pazos, and D. D. Sampson, “Review of tissue simulating phantoms with controllable optical, mechanical and structural properties for use in optical coherence tomography,” Biomed. Opt. Express 3(6), 1381–1398 (2012). [CrossRef]  

39. Q. Wu and A. T. Yeh, “Rabbit cornea microstructure response to changes in intraocular pressure visualized by using nonlinear optical microscopy,” Cornea 27(2), 202–208 (2008). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (6)

Fig. 1.
Fig. 1. Illustrative OCT scans from each considered study with ROI marked with red lines. Cyan lines indicate phantom borders (top image) or epithelium, Bowman’s layer, and endothelium (the middle and bottom image).
Fig. 2.
Fig. 2. Distance values for phantoms with decreasing concentration of blue dye particles. Illustrative images inserted at the bottom of the plot present the ROIs of the OCT images of the phantoms for the concentrations C8, C5, and C2.
Fig. 3.
Fig. 3. Estimated distances as functions of ROI area for the phantom of concentration C9. S denotes the original ROI area of $600\times 220$ pixels.
Fig. 4.
Fig. 4. Values of parameters of the considered distributions for decreasing concentrations of dye particles in the phantom study.
Fig. 5.
Fig. 5. Boxplots of $D_{\mathrm {MSE}}$, $D_{\mathrm {KS}}$, $D_{\mathrm {MMD}}$, and $D_{\mathrm {CR}}$ values for Experiment 1 (examining the impact of IOP) and Experiment 2 (examining the impact of its duration) on porcine corneas (ex-vivo study). Asterisks indicate statistical significance of differences between the values of statistical distances for adjacent IOP levels (blue ones) or time points (red ones), assessed using paired t-test (* $p<0.05$; ** $p<0.01$; *** $p<0.001$).
Fig. 6.
Fig. 6. Means and standard errors of the goodness-of-fit (GoF) for a set of different distributions for the OCT images of porcine eyeballs in the ex-vivo study in Experiment 1.

Tables (1)

Tables Icon

Table 1. The results of repeated measures ANOVA for two ex-vivo porcine eye experiments and the results of correlation analysis for the in-vivo human eye study for both non-parametric and parametric approaches.

Equations (14)

Equations on this page are rendered with MathJax. Learn more.

p A ( A ) = A σ 2 e A 2 2 σ 2 , A 0.
σ ^ M L E = 1 2 n i = 1 n A i 2 ,
σ ^ M L E = 1 2 n 1 1 n i = 1 n A i 2 i = 1 n A i 2 = 2 2 ,
σ ^ B a y e s = 2 Γ ( n + 2 ) 2 Γ ( n + 5 2 ) i = 1 n A i 2 .
σ ^ B a y e s = 2 Γ ( n + 2 ) 2 Γ ( n + 5 2 ) 1 1 n i = 1 n A i 2 i = 1 n A i 2 = 2 2 Γ ( n + 2 ) Γ ( n + 5 2 ) n ,
p A ( A ) = 4 Γ ( α ) α A 2 ( α A 2 A 2 ) α 2 K α 1 ( 2 α A 2 A 2 ) ,
e C D F = 1 n i = 1 n 1 { X i t } ,
K D E = 1 n i = 1 n 1 h K ( x X i h ) ,
e C F = 1 n i = 1 n e j t X i , j = 1 .
C R = σ ^ X X ¯ ,
D K S = sup x | e C D F ( x ) C D F B D ( x ) |
D M S E = 1 n i = 1 n ( K D E ( x i ) P D F B D ( x i ) ) 2
D M M D = 1 n i = 1 n e C F ( x i ) 1 n i = 1 n C F B D ( x i )
D C R = C R 0.5227
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.