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

Speckle statistics of biological tissues in optical coherence tomography

Open Access Open Access

Abstract

The speckle statistics of optical coherence tomography images of biological tissue have been studied using several historical probability density functions. Here, we propose a new theoretical framework based on power-law functions, where we hypothesize that an underlying power-law distribution governs scattering from tissues. Thus, multi-scale scattering sites including the fractal branching vasculature will contribute to power-law probability distributions of speckle statistics. Specifically, these are the Burr type XII distribution for speckle amplitude, the Lomax distribution for intensity, and the generalized logistic distribution for log amplitude. Experimentally, these three distributions are fitted to histogram data from nine optical coherence tomography scans of various samples and biological tissues, in vivo and ex vivo. The distributions are also compared with classical models such as the Rayleigh, K, and gamma distributions. The results indicate that across OCT datasets of various tissue types, the proposed power-law distributions are more appropriate models yielding novel parameters for characterizing the physics of scattering from biological tissue. Thus, the overall framework brings to the field new biomarkers from OCT measures of speckle in tissues, grounded in basic biophysics and with wide applications to diagnostic imaging in clinical use.

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

1. Introduction

Speckle is a granular pattern seen in signals or images that is caused by the interference of coherent waves with random phases and known or random amplitudes. [1] The study of speckle phenomena has a long history dating back to the time of Isaac Newton, with an increasing multitude of recent applications in optics, radar, and ultrasound. Optical coherence tomography (OCT) and ultrasound are two medical imaging modalities with prominent speckle. For some applications such as high-resolution imaging, speckle is considered to be undesirable noise and many studies attempt to eliminate its presence. [210] Other studies choose to utilize speckle for physical modeling or characterization of tissue samples. [1116] Studies of speckle amplitude statistics in acoustics and optics have led to the usage of various probability density functions (PDFs) such as the Rayleigh distribution, the K distribution, the Rice distribution, gamma distributions, and many others. [1722]

Although OCT is an interferometric technique and ultrasound utilizes time-of-flight measurements, the mathematics describing wave propagation and wave phenomena such as speckle can be applicable to both acoustical and optical imaging modalities. Recently, a new model hypothesized that in normal soft tissue, the dominant scattering elements are cylinders from fractal branching vasculature. [2326] The results apply generally to multi-scale distributions of scattererers and three new and distinct probability distributions were found to characterize ultrasound speckle in biological tissue. These were the Burr type XII distribution, the Lomax distribution, and the generalized logistic distribution.

In this paper, this novel framework for speckle is extended to OCT scans of various biological tissues. Metrics for assessing appropriate regions of interest (ROIs) and evaluating the statistical validity of the distributions are also presented. A comparison of the new distributions with distributions found in the literature are presented based on clinically relevant ex vivo and in vivo samples, and the power-law distributions are shown to provide a superior description of the speckle statistics. Finally, the key power-law parameter estimated from speckle is capable of discriminating between tissues and thus introduces a new biomarker for clinical use.

2. Theory

2.1 Modeling of biological tissue

Parker et al. [26] recently derived the first order speckle statistics of biological tissue in ultrasound imaging under the assumptions of weak scattering (using the Born approximation) originating from fractal branching of vasculature represented by cylinders. This derivation leads to power-law functions that dictate the PDFs for the echo amplitude and intensity histograms. One of the properties about power laws is their invariance to scale, which allows us to bridge between ultrasound and OCT, two modalities of different resolutions. The framework for the derivation is also directly applicable to multi-scale optical scattering within OCT and is summarized here.

First, consider a distribution of scattering structures, from large to small, within the volume. We assume the distribution follows a power-law distribution in size, with relatively fewer larger scatterers. A power-law distribution and spatial correlation function are also consistent with generalized fractal models. [27] Thus, in scanning a volume, the probability of encountering a scatterer of characteristic dimension $a$ is given as [26,28]

$$p(a) = \frac{b-1}{a_{\textrm{min}}} \left(\frac{a}{a_{\textrm{min}}}\right)^{{-}b}$$
where $b$ is the power-law coefficient representing the multi-scale nature of the tissue structures, and $a_{\textrm {min}}$ represents the minimum size or lower limit of dimensions of the scattering structures that are detectable. Previous studies have shown that variations in the index of refraction within tissues obey a power law down to the sub-micron scale. [29] Thus, this model is appropriate for tissues.

Secondly, we assume each scatterer of dimension $a$ produces a detected amplitude $A$ or intensity $I$ according to the theory of backscattered waves. In general, canonical scattering elements such as spheres and cylinders have been characterized by power series solutions. [3032] The dependence of backscatter on frequency and dimension is complicated, but can be characterized by well-known long-wavelength, short-wavelength, and transition or Mie scattering regimes. [33] However, in the sub-resolvable regime where the scatterer is smaller than a wavelength, backscatter models exhibit similar trends as scatterer size increases from the lower limit (long wavelength or Rayleigh backscatter regime) towards larger sizes approaching the wavelength. Across this sub-resolvable range, the scattering versus size curves increase progressively until the Mie scattering regime is reached. As a first-order approximation to this behavior, we have employed a linear offset model covering the sub-resolvable regime where

$$I(a) = I_0 (a-a_{\textrm{min}})$$
where both $I_0$ and $a_{\textrm {min}}$ are dependent on system parameters such as wavelength and gain, and the lower limit of system detectability, including quantization and noise floor. With this linear monotonic function, the probability of occurrence can be simply mapped into the probability of amplitude or intensity using the probability transformation rule. [24,26]

2.2 Probability distributions for amplitude, intensity, and log of amplitude

The PDF for the histogram of amplitudes $x$ is derived by Parker et al. [26] and is given by

$$p(x;d,b) = \frac{2x(b-1)}{d^2 \left[ \left( \frac{x}{d} \right)^2 + 1 \right]^b}, \quad x>0$$
which is a special case of the Burr type XII distribution as a two-parameter fit. The parameters are the power-law exponent $b$ and a scale factor $d$. When using this PDF for fitting OCT speckle amplitude, it is convenient to normalize by setting $x=\frac {A}{\sqrt { \langle A^2 \rangle }}$, where the denominator represents the root mean square (RMS) value. [18,21] Otherwise, the normalization constant can be incorporated into the $d$ parameter.

The PDF for the histogram of intensities is given by

$$p(x;d,b) = \frac{(b-1)d^{b-1}}{(x+d)^b}, \quad x>0$$
which is the Lomax distribution or Pareto type II distribution. When using this PDF for fitting OCT speckle intensity $I=|A|^2$, it is convenient to normalize by setting $x=\frac {I}{\langle I \rangle }$, where $\langle I \rangle$ is the mean intensity.

Finally, the PDF for the histogram of log amplitude defined by $y=\log (A)$ is given by

$$p(y;d,b) = \frac{2(b-1)\exp{(2y)}}{d^2 \left[ \frac{\exp{(2y)}}{d^2} + 1 \right]^b}, \quad -\infty<y<\infty$$
which is a transformed version of the generalized logistic type I distribution. Typically, in OCT, it is beneficial to display an image of the log of the amplitude or intensity to better visualize the dynamic range. Thus, this distribution is useful in capturing the histograms of most conventional display values. Furthermore, these three PDFs are all well characterized in the statistics and econometrics literature, with known cumulative distribution functions and moments. [34]

2.3 Theoretical importance of the exponent parameter

The above three PDFs all contain a power-law or exponent parameter $b$. In power-law and related functions, the exponent parameter is a valuable parameter of interest in many applications. [28] In this paper’s context, the exponent parameter $b$ is important in tissue characterization.

According to Carroll-Nellenback, et al., [35] a simple fractal distribution of vessels within normal tissues would provide a value of approximately $b=2.7$. However, the number of scatterers per sample volume and the index of refraction can increase the exponent parameter. Thus, the exponent parameter provides information about the tissue’s scattering properties and structure, and acts as a biomarker for differentiating different tissues.

2.4 Historical probability distributions for comparison

In the literature, there are other PDFs used to model OCT speckle amplitudes and intensities based on consideration of random point scatterers or more complex distributions from radar and other fields. The most prevalent of these is the Rayleigh distribution for speckle amplitude and the exponential distribution for speckle intensity, which are given by: [18,20]

$$\begin{aligned}p(A) &= \frac{2A}{\sqrt{\langle A^2 \rangle}} \exp{\left( -\frac{A^2}{\langle A^2 \rangle} \right)}, \quad A>0 \\ p(I) &= \frac{1}{\langle I \rangle} \exp{\left( -\frac{I}{\langle I \rangle} \right)}, \quad I>0 \end{aligned}$$
The Rayleigh distribution is suitable for the case of a large number of scatterers in a homogeneous medium.

Another distribution used by Weatherbee et al. is the K distribution, which has also been explored in ultrasound and radar. [18,21] The K distribution is modeled for the case of a small number of scatterers, and the PDFs for the amplitudes and intensities are given by

$$\begin{aligned}p(A;\alpha) &= \frac{4}{\Gamma(\alpha)} \sqrt{\frac{\alpha}{\langle A^2 \rangle}} \left( \frac{\alpha A^2}{\langle A^2 \rangle} \right)^{\alpha/2} K_{\alpha-1}\left( 2 \sqrt{\frac{\alpha A^2}{\langle A^2 \rangle}} \right), \quad A>0 \\ p(I;\alpha) &= \frac{2}{\Gamma(\alpha)} \sqrt{\frac{\alpha}{I \langle I \rangle}} \left( \frac{\alpha I}{\langle I \rangle} \right)^{\alpha/2} K_{\alpha-1}\left( 2 \sqrt{\frac{\alpha I}{\langle I \rangle}} \right), \quad I>0 \end{aligned}$$
where $\Gamma (\cdot )$ is the gamma function and $K(\cdot )$ is a modified Bessel function of the second kind, and $\alpha$ is the shape parameter.

A third distribution for comparison is the gamma distribution, as used by Kirillin et al. for modeling speckle amplitude, and is given by [16]

$$\begin{aligned}p(A;\alpha,\beta) &= \frac{1}{\Gamma(\alpha)} \beta^\alpha A^{\alpha-1} \exp{\left( -\beta A \right)}, \quad A>0 \\ p(I;\alpha,\beta) &= \frac{1}{2\Gamma(\alpha)} \beta^\alpha I^{\frac{\alpha}{2}-1} \exp{\left( -\beta \sqrt{I} \right)}, \quad I>0 \end{aligned}$$
where $\alpha$ and $\beta$ are two shape parameters.

3. Methods

3.1 OCT scans of various biological tissues

A swept source OCT (SS-OCT) system is used to scan various biological tissue. It is implemented with a swept source laser (HSL-2100-WR, Santec, Aichi, Japan) with a center wavelength of 1310 nm and full-width half-maximum (FWHM) bandwidth of 170 nm. The lateral resolution is 20 µm and the FWHM of the axial point spread function after dispersion compensation is 8 µm in air. The maximum sensitivity of the system was measured to be 120 dB. The imaging depth was measured to be 5 mm in air (-10 dB sensitivity roll-off). The scanning lens has an $f$-number of $f/9$. The SS-OCT system is controlled with LabVIEW (Version 14, National Instruments, Austin, Texas, USA).

The following tissues were scanned and analyzed with the SS-OCT system: mouse brain and liver, pig brain and cornea, and chicken muscle all ex vivo as well as human hand (skin) in vivo. In addition, two gelatin phantoms (5% with and without milk for optical scattering) were also scanned and analyzed. Milk in the form of coffee creamer was added as 0.2% of phantom weight. Additional scans of three pig corneas were performed to assess reproducibility. The number of A-lines for each scan was either 100, 500, or 1000. Variations in the number of A-lines do not change the speckle statistics, as long as an adequate number of pixels are used (e.g. greater than 1,000 pixels, which is easily covered by a 10 $\times$ 100 ROI) over an appropriate field of view (e.g. 5-10 mm). A-lines were not averaged for analysis to preserve the statistics of power-law related tail behaviors. Focal planes were set at approximately 0.2 mm below the surface of the sample apex. Amplitude, intensity, and log amplitude histograms are generated from specific ROIs in these samples.

3.2 Fitting distributions using maximum likelihood estimation

The distributions specified in Sections 2.2 and 2.4 are fitted to respective amplitude and intensity histograms using maximum likelihood estimation (MLE). This technique estimates the parameters of a PDF by maximizing a known likelihood or log-likelihood function. [36] The iterative maximization algorithm for MLE and all other analysis aspects were conducted in MATLAB 2020b (MathWorks, Natick, Massachusetts, USA). Curve fitting using an alternative least-squares approach or related methods can result in systemic errors, especially when estimating a power-law or exponent parameter. [37] This is due to a combination of factors such as violations of the underlying assumptions when using these curve fitting methods and variability from histogram binning methods. Therefore, MLE is a more accurate approach to this paper’s studies.

3.3 Metric for specifying an appropriate region of interest

Relative uniformity across the ROI is an important, yet difficult to quantify, requirement for assessing appropriate speckle statistics. Attenuation along depth and shadowing effects are two examples of phenomena that would reduce the validity of an ROI for appropriate speckle statistics. Within a region of macroscopically similar tissue, we assume spatial stationarity and ergodicity of speckle statistics, and therefore the average along any projection line would have an identical expected value. This forms the basis of our following ad hoc metric for accepting or rejecting an ROI for analysis.

Simulations were performed to obtain a threshold for selecting an ROI with acceptably uniform statistics. In Fig. 1(a), a simulated OCT B-mode image (500 $\times$ 500 pixels) using a Burr type XII distribution with parameters $d=10$ and $b=3$ is shown with attenuation and shadowing effects. Figure 1(b) shows how spatial integration profiles vary along the axial and lateral directions, which can be visually correlated to the B-mode image. Specifically, integrations of speckle amplitude are performed along the two directions using intervals of 5 pixels and are normalized so that the two profiles can be displayed on the same scale. Using a Monte Carlo method, 10,000 ROIs of the B-mode image are randomly generated with a minimal size of 15 $\times$ 15 pixels to ensure adequate statistics. The resulting speckle statistics were fitted to a Burr type XII distribution using MLE. The ROIs were quantified by the maximum percent change (i.e. ratio of maximum difference over minimum value) in their spatial integration profiles, $\Delta$. Figure 1(c) shows how the estimated parameter $\hat {b}$ varies as a function of the ROI’s $\Delta$. Using this simulation, we established an ad hoc requirement that any appropriate ROI’s $\Delta$ should not exceed a 50% change in their spatial integration profiles laterally and axially.

 figure: Fig. 1.

Fig. 1. Determining an ROI metric. (a) Simulated B-mode image using a Burr type XII distribution ($d=10,\,b=3$) with attenuation and shadowing. (b) Depth and lateral profiles based on spatial integration laterally and along depth, respectively. (c) Multiple ROIs and MLE fitting were used to calculate $\hat {b}$ as a function of maximum percent change in profiles, $\Delta$.

Download Full Size | PDF

In addition, the estimated attenuation coefficient, denoted as $\hat {\mu }$, is estimated using a depth-resolved approach, as demonstrated by Vermeer, et al. [38]

$$\hat{\mu} = \frac{1}{2\,dz} \log{\left( 1+\frac{I[i]}{\sum_{i+1}^{\infty} I[i]} \right)}$$
where $i$ denotes the $i^{\textrm {th}}$ pixel, $I$ is the intensity, and $dz$ is the pixel size along depth. The average attenuation coefficient, denoted as $\hat {\mu }_{\textrm {avg}}$, along with the standard deviation are computed for each sample ROI. The attenuation coefficient provides an additional parameter for tissue characterization. The Pearson correlation coefficient $r$ between the attenuation coefficient and the exponent parameter is also reported as $r(\textrm {degrees of freedom}) = r \textrm { statistic}, p \textrm { value}.$

3.4 Comparing multiple distributions

There are numerous methods for direct comparison of models. These methods include the Kolmogorov-Smirnov (KS) test, the Anderson-Darling (AD) test, the likelihood ratio test (LRT), and the Akaike information criterion (AIC). [37,39,40] The KS and AD test statistics provide $p$-values and indicate which models are best rejected while simultaneously indicating which model is the best fit. The AD test is similar to the KS test but gives more weight to distribution tails. However, it can be too conservative with estimates of power-law functions and requires customized tests for the PDFs described in Sections 2.2 and 2.4. Instead, we follow the rapid KS procedure as outlined in Clauset et al. [37] By performing the KS test on the sample data and simulated sample distributions from the MLE fit, a $p$-value can be estimated by using a Monte Carlo method. In this case, a very large $p$-value (> 0.9) indicates a good fit, while a small value indicates that the model is not an appropriate one. When comparing models, the one with the largest KS $p$-value can be considered the best, although the best model may still not be a good fit.

The second test used is the LRT. The LRT generates a ratio of the log likelihoods of two models for comparison, which we denote by $R$. [37,39] In this case, the LRT is used to compare the distributions in Section 2.4 (Rayleigh/Exponential, K, and Gamma distributions) with the Burr/Lomax distributions. The sign of $R$ indicates which distribution is the better fit: if $R<0$, then the Burr/Lomax distribution is the preferred distribution. A $p$-value is also estimated with $R$ to acknowledge that a true $R=0$ may fluctuate either positively or negatively. If the $p$-value is small, then the sign of $R$ is a good indicator of which distribution is the better fit. Otherwise for a non-small $p$-value (> 0.1), the LRT is inconclusive.

The final metric used is the AIC, which is defined as

$$\textrm{AIC} = 2k - 2\log\left( \hat{L} \right)$$
where $k$ is the number of estimated parameters and $\hat {L}$ is the model’s likelihood. The AIC provides a relative measure of a model when compared to other models, where a smaller value indicates a better model. [40] These three methods (KS test, LRT, and AIC) are used to determine which distribution fits best for speckle amplitude (Burr type XII, Rayleigh, K, or Gamma) and intensity (Lomax, Exponential, K, or Gamma).

4. Results

4.1 Evaluating the Burr type XII, Lomax, and generalized logistic distributions

Figure 2 demonstrates the speckle analysis performed for a mouse brain sample. Figure 2(a) shows the B-mode image and the ROI used for the speckle statistics, and Fig. 2(b) verifies that this ROI is appropriate using the metric defined in Section 3.3. The MLE fits for the amplitude values using the Burr type XII distribution are shown in Fig. 2(c)-(d), on both a linear and logarithmic scale. Figure 2(e)-(f) reports the MLE fits for the intensity values using the Lomax distribution, and Fig. 2(g)-(h) shows the MLE fits for the log of amplitude values using the generalized logistic distribution. All MLE fits to the prescribed PDFs result in a consistent exponent parameter estimate $\hat {b}=6.06$ (95% Confidence Interval: $[5.83,\,6.28]$) for mouse brain. Figure 3 demonstrates the same type of analysis performed for mouse liver, which results in a consistent estimate of $\hat {b}=5.95$ (95% Confidence Interval: $[5.70,\,6.21]$). Visual inspection indicates that the Burr type XII, the Lomax, and the generalized logistic distributions are good fits to the data, but statistical analyses of these distributions are further investigated in the next section.

 figure: Fig. 2.

Fig. 2. Mouse brain with cranial window. (a) Single unfiltered and unaveraged B-mode frame with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles demonstrating that the maximum percent change $\Delta$ does not exceed 50%. (c) MLE fit to the Burr type XII distribution in a linear scale with linear binning of amplitude histogram data. (d) Same as (c) but in a log-log scale with logarithmic binning for visualization of tail behavior. (e) MLE fit to the Lomax distribution in a linear scale with linear binning of intensity histogram data. (f) Same as (d) but in a log-log scale with logarithmic binning to visualize tail behavior. (g) MLE fit to the generalized logistic distribution in a linear scale with linear binning of log amplitude histogram data. (h) Same as (g) but in log-log scale with logarithmic binning. All MLE fits result in an exponent parameter of $\hat {b}=6.06$ (95% Confidence Interval: $[5.83,\,6.28]$) for mouse brain.

Download Full Size | PDF

 figure: Fig. 3.

Fig. 3. Excised mouse liver in phosphate-buffered saline. (a) Single unfiltered and unaveraged B-mode frame with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles demonstrating that the maximum percent change $\Delta$ does not exceed 50%. (c) MLE fit to the Burr type XII distribution in a linear scale with linear binning of amplitude histogram data. (d) Same as (c) but in a log-log scale with logarithmic binning for visualization of tail behavior. (e) MLE fit to the Lomax distribution in a linear scale with linear binning of intensity histogram data. (f) Same as (d) but in a log-log scale with logarithmic binning to visualize tail behavior. (g) MLE fit to the generalized logistic distribution in a linear scale with linear binning of log amplitude histogram data. (h) Same as (g) but in log-log scale with logarithmic binning. All MLE fits result in an exponent parameter of $\hat {b}=5.95$ (95% Confidence Interval: $[5.70,\,6.21]$) for mouse liver.

Download Full Size | PDF

4.2 Comparing multiple distributions

Figure 4 demonstrates multiple distribution fits for a pig cornea sample. Figure 4(a) shows the B-mode frame for the pig cornea with an appropriate ROI as quantified in Fig. 4(b). Figure 4(c)-(d) depicts the MLE fits for amplitude histogram data using the Burr type XII, the Rayleigh, the K, and the gamma distributions as specified in Section 2.4 in both a linear and logarithmic scale. The intensity histogram data and the Lomax, exponential, K, and gamma distribution MLE fits are shown in Fig. 4(e)-(f). All respective MLE fits for amplitude and intensity result in consistent parameter estimates (e.g. $\hat {\alpha }$ and $\hat {\beta }$ are the same for the gamma distributions for both amplitude and intensity). Figure 5 shows the sample analysis done for a human hand’s back. Once again, the MLE parameter estimates are consistent between amplitude and intensity fits for each pair of distributions (Burr/Lomax, Rayleigh/exponential, K, or gamma). The complete set of data for all tissue or phantom samples taken can be found in Supplemental Fig. S1-S9. In some cases, on the logarithmic scale, there are slight deviations between the data points and the estimated Burr and Lomax distributions. Since there are a finite number of pixels in a single ROI, there is a minimum probability estimate that can be captured by histogram data. Additionally, the bins in the tail carry a small number of samples, and so statistical fluctuations account for a larger proportion of these samples. Hence, results in the tail become noisy, and deviations from the estimated distributions are expected.

 figure: Fig. 4.

Fig. 4. Pig cornea. (a) Single unfiltered and unaveraged B-mode image with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles that do not exceed a $\Delta$ of 50%. (c) MLE fits to normalized amplitude with the Burr, Rayleigh, K, and Gamma distributions on a linear scale with linear binning. (d) Same as part (c) except on a log-log scale with logarithmic binning to visualize tail behavior. (e) MLE fits to normalized intensity with the Lomax, Exponential, K, and Gamma distributions on a linear scale with linear binning. (f) Same as part (e) except on a log-log scale with logarithmic binning to visualize tail behavior.

Download Full Size | PDF

 figure: Fig. 5.

Fig. 5. Human hand (back side). (a) Single unfiltered and unaveraged B-mode image with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles that do not exceed a $\Delta$ of 50%. (c) MLE fits to normalized amplitude with the Burr, Rayleigh, K, and Gamma distributions on a linear scale with linear binning. (d) Same as part (c) except on a log-log scale with logarithmic binning to visualize tail behavior. (e) MLE fits to normalized intensity with the Lomax, Exponential, K, and Gamma distributions on a linear scale with linear binning. (f) Same as part (e) except on a log-log scale with logarithmic binning to visualize tail behavior.

Download Full Size | PDF

The next step is to perform the statistical tests and metrics for comparing multiple distributions described in Section 3.4. The results of all data samples are detailed in Supplemental Tables S1-S9, but a summary table of the comparison outcomes is provided in Supplemental Table S10. In all samples except for one, the best model for amplitude was the Burr type XII distribution and the best model for intensity was the Lomax distribution. The one exception is for pig brain, in which the K distribution was the best fit for both amplitude and intensity data.

4.3 Estimating the exponent parameter of multiple types of biological tissues

The Burr type XII, Lomax, and generalized logistic distributions all have two parameters: $d$ and $b$. $d$ can be seen as a system-dependent (and gain-dependent) normalization parameter, but the exponent parameter $b$ is of interest for tissue characterization. For all samples, $\hat {b}$ and its 95% confidence interval were obtained via MLE in Supplemental Fig. S1-S9. A summary of the exponent values for various samples types are shown in Table 1, in order of increasing $\hat {b}$, along with average attenuation coefficients. Note that for all samples except the pig brain, the Burr/Lomax distribution was determined to be the best model in the previous section. Reproducibility of the method is demonstrated in Supplemental Figures S10-S12 and Supplemental Tables S11-S14. Using the three extra pig corneas and the one presented above, the average and standard deviation of the exponent parameter was $3.00 \pm 0.04$ with a combined 95% confidence interval range of $[2.92,\,3.10]$.

Tables Icon

Table 1. Estimated exponent parameter values along with the 95% confidence interval for all samples. The average attenuation coefficient $\hat {\mu }_{\textrm {avg}}$ in the ROI is also reported with the standard deviation. CI - confidence interval.

5. Discussion and conclusion

In Section 4.1, MLE fits to the Burr type XII, Lomax, and generalized logistic distributions show that these three PDFs fit well to the amplitude, intensity, and log amplitude histogram data for two sample tissues (mouse brain and liver). The estimated exponential parameter $\hat {b}$ is also universally consistent due to the framework of MLE. Thus, any one of the three PDFs can be reasonably used to estimate the exponential parameter $b$. In Section 4.2, we further demonstrate the merits of using these PDFs (Burr/Lomax), as they are statistically compared with conventional PDFs described in the literature (Rayleigh/exponential, K, and gamma). Out of nine total samples, eight of them demonstrated that Burr/Lomax are indeed the best fits for the amplitude and intensity histogram data. In the one exception (pig brain), the differences between the Burr/Lomax fits and the K fits were statistically significant but relatively small. Another trend that was noticed was the fact that if the Burr/Lomax estimated a relatively high $\hat {b} \sim 6$, then the K and gamma fits were relatively closer to the Burr/Lomax fits, and they are reasonable fits to the data (this can be verified by checking that the KS $p$-values are large). On the other hand, if the Burr/Lomax’s exponent parameter was small, e.g. $\hat {b} \sim 3$, then the K and gamma distributions were poorer fits to the data (small KS $p$-values). Overall, the Burr/Lomax distributions indicate that they are reliable and generalizable to characterizing the speckle statistics of various biological tissue in OCT.

In Section 4.3, the exponent parameter $\hat {b}$ is tabulated and compared across multiple biological tissues. The range for $\hat {b}$ is from $3 \sim 6$. There may be some reasonable explanations for the general trend of these values. The samples with lower values of $\hat {b} \sim 3$ are all tissues that are relatively transparent to light (eye, skin, gelatin), whereas the samples with higher values of $\hat {b} \sim 6$ are considered optically denser tissues (liver, brain). The differences between the phantom studies also supports this concept. The gelatin phantom with milk (an added optical scatterer) has a higher $\hat {b}$ value. Furthermore, since the Burr/Lomax distributions are much better fits to the speckle statistics than Rayleigh/K/gamma distributions for tissues with low values of $\hat {b}$, then the Burr/Lomax may be even more important in characterizing these types of tissues. OCT currently has strong clinical applications in ophthalmology and dermatology, and so a framework involving $\hat {b}$ would provide a useful biomarker for pathology or disease.

In addition, there is a moderate positive correlation between the estimated exponent parameter $\hat {b}$ and the average attenuation coefficient $\hat {\mu }_{\textrm {avg}}$, $r(7) = 0.52, p = 0.274$. The high $p$-value is due to the small sample size. The reported attenuation coefficients are consistent with some other values reported in the literature. [38,4146] Attenuation is contributed by both scattering and absorption, and thus one would expect that optically denser tissues with high absorption and/or scattering would lead to an increased exponent parameter $b$ as well as a higher attenuation coefficient $\mu$. This can be seen in Table 1 generally, with the human hand (back) sample as a notable outlier. Furthermore, future studies would be needed to correlate the speckle parameter $b$ with other common scattering parameters such as the mean free path and transport mean free path. Further insight on these new PDFs may also allow for characterization of wave transport and light scattering in the context of disordered systems.

Our test for spatial uniformity does not guarantee stationary statistics within the ROI, and so we have also examined the distribution of local second moments as proposed by Bromberg and Cao, in their supplemental material. [47] In this examination, the local contrast is defined as $C=\sqrt {\langle I^2 \rangle / \langle I \rangle ^2 - 1}$ which would be everywhere close to unity for a stationary Rayleigh process. We found that in our samples of mouse brain, mouse liver, pig cornea, and human hand (back), the local contrast $C$ as a function of the number of data points or speckles in the ROI converges to the following values: 1.28, 1.29, 2.00, and 1.84, respectively. These were also found to have small (less than 10%) local variations within the ROI for sub-regions (with at least 10 speckle areas) and together these results imply stationarity, and that the processes are super-Rayleigh.

Another potential concern that has been investigated is the effect of correcting the image for the confocal point spread function and sensitivity roll-off function. When applying the estimates for these corrections using the methods demonstrated by Fiske, et al., [46] it was found that these system corrections typically alter the Burr/Lomax exponent parameters by less than 0.01 and the attenuation coefficients by less than 0.1. The magnitudes of these changes are well within the 95% confidence intervals of the exponent parameters and the standard deviations of the average attenuation coefficients for all samples. Furthermore, we have shown that the ad hoc test ensures that depth-dependent effects of attenuation are limited, and also that local statistics are not spatially varying. Thus, the estimates of attenuation within the selected ROIs are dominated by tissue attenuation, as opposed to the confocal point spread function and the sensitivity roll-off function. In addition, these corrections are only approximate and not exact. However, detailed analysis of these effects on estimating the power-law parameter $b$ in larger ROIs should be conducted in future research.

In summary, this study has presented three new PDFs for fitting speckle statistics in various biological tissues: the Burr type XII distribution for amplitude, the Lomax distribution for intensity, and the generalized logistic distribution for log of amplitude. Furthermore, an ad hoc metric for verifying an appropriate ROI, the MLE fitting technique, and methodology for statistical comparison of multiple distributions are also presented. Future work is needed to closely link the speckle parameter $b$ with spatial measures of tissue microstructure in 3D. [29,35,48] In addition, the $b$ parameter relates to power-law models of scattering from tissues, which have been explored previously, but future studies are needed to relate $b$ to other common scattering parameters. [49,50] These independent measures on a larger set of tissues will be needed to verify that the results of this study are reproducible across tissues, and to determine the importance of the exponent parameter $b$ in multiple biological and clinical applications.

Funding

National Institutes of Health (F30AG069293, R21AG070331, R21EB025290).

Acknowledgments

Gary Ge is supported by the National Institute on Aging of the National Institutes of Health under award number F30AG069293. The work is also supported by NIH grants R21EB025290 and R21AG070331. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors would also like to thank Fernando Zvietcovich for providing the OCT scans of pig brain and chicken muscle.

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. Correspondence should be addressed to the corresponding author.

Supplemental document

See Supplement 1 for supporting content.

References

1. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (SPIE The International Society for Optics and Photonics, 2020), 2nd ed.

2. M. Bashkansky and J. Reintjes, “Statistics and reduction of speckle in optical coherence tomography,” Opt. Lett. 25(8), 545–547 (2000). [CrossRef]  

3. K. Z. Abd-Elmoniem, A. M. Youssef, and Y. M. Kadah, “Real-time speckle reduction and coherence enhancement in ultrasound imaging via nonlinear anisotropic diffusion,” IEEE Trans. Biomed. Eng. 49(9), 997–1014 (2002). [CrossRef]  

4. P. Michael, G. Erich, A. L. Rainer, F. Adolf Friedrich, and K. H. Christoph, “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt. 8(3), 565–569 (2003). [CrossRef]  

5. D. C. Adler, T. H. Ko, and J. G. Fujimoto, “Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter,” Opt. Lett. 29(24), 2878–2880 (2004). [CrossRef]  

6. J. Kim, D. Miller, E. Kim, S. Oh, J. H. Oh, and T. Milner, “Optical coherence tomography speckle reduction by a partially spatially coherent source,” J. Biomed. Opt. 10(6), 064034 (2005). [CrossRef]  

7. A. Ozcan, A. Bilenca, A. E. Desjardins, B. E. Bouma, and G. J. Tearney, “Speckle reduction in optical coherence tomography images using digital filtering,” J. Opt. Soc. Am. A 24(7), 1901–1910 (2007). [CrossRef]  

8. Y. Guo, H. D. Cheng, J. Tian, and Y. Zhang, “A novel approach to speckle reduction in ultrasound imaging,” Ultrasound Medicine & Biol. 35(4), 628–640 (2009). [CrossRef]  

9. J. Park, J. B. Kang, J. H. Chang, and Y. Yoo, “Speckle reduction techniques in medical ultrasound imaging,” Biomed. Eng. Lett. 4(1), 32–40 (2014). [CrossRef]  

10. O. Liba, M. D. Lew, E. D. SoRelle, R. Dutta, D. Sen, D. M. Moshfeghi, S. Chu, and A. de la Zerda, “Speckle-modulating optical coherence tomography in living mice and humans,” Nat. Commun. 8(1), 15845 (2017). [CrossRef]  

11. T. A. Tuthill, R. H. Sperry, and K. J. Parker, “Deviations from Rayleigh statistics in ultrasonic speckle,” Ultrason. Imaging 10(2), 81–89 (1988). [CrossRef]  

12. J. M. Thijssen, “Ultrasonic speckle formation, analysis and processing applied to tissue characterization,” Pattern Recognit. Lett. 24(4-5), 659–675 (2003). [CrossRef]  

13. B. Karamata, K. Hassler, M. Laubscher, and T. Lasser, “Speckle statistics in optical coherence tomography,” J. Opt. Soc. Am. A 22(4), 593–596 (2005). [CrossRef]  

14. L. Di, N. Rao, C. Kuo, S. Bhatt, and V. Dogra, “Independent component analysis applied to ultrasound speckle texture analysis and tissue characterization,” in 2007 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, (2007), pp. 6523–6526.

15. 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]  

16. 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]  

17. O. S. Al-Kadi, D. Y. F. Chung, C. C. Coussios, and J. A. Noble, “Heterogeneous tissue characterization using ultrasound: A comparison of fractal analysis backscatter models on liver tumors,” Ultrasound Medicine & Biol. 42(7), 1612–1626 (2016). [CrossRef]  

18. 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–2730 (2016). [CrossRef]  

19. J. Kalkman, “Fourier-domain optical coherence tomography signal analysis and numerical modeling,” Int. J. Opt. 2017, 1–16 (2017). [CrossRef]  

20. 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]  

21. 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]  

22. V. Y. Zaitsev, L. A. Matveev, A. L. Matveyev, G. V. Gelikonov, and V. M. 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]  

23. K. J. Parker, “Shapes and distributions of soft tissue scatterers,” Phys. Med. Biol. 64(17), 175022 (2019). [CrossRef]  

24. K. J. Parker, “The first order statistics of backscatter from the fractal branching vasculature,” J. Acoust. Soc. Am. 146(5), 3318–3326 (2019). [CrossRef]  

25. K. J. Parker and S. S. Poul, “Speckle from branching vasculature: dependence on number density,” J. Med. Imaging 7(2), 027001 (2020). [CrossRef]  

26. K. J. Parker and S. S. Poul, “Burr, Lomax, Pareto, and logistic distributions from ultrasound speckle,” Ultrason. Imaging 42(4-5), 203–212 (2020). [CrossRef]  

27. T. Vicsek, Fractal Growth Phenomena (World Scientific, 1992).

28. M. E. J. Newman, “Power laws, pareto distributions and Zipf’s law,” Contemp. Phys. 46(5), 323–351 (2005). [CrossRef]  

29. J. M. Schmitt and G. Kumar, “Turbulent nature of refractive-index variations in biological tissue,” Opt. Lett. 21(16), 1310–1312 (1996). [CrossRef]  

30. L. Rayleigh, “X. on the electromagnetic theory of light,” The London, Edinburgh, Dublin Philos. Mag. J. Sci. 12(73), 81–101 (1881). [CrossRef]  

31. G. Mie, “Beiträge zur optik trüber medien, speziell kolloidaler metallösungen,” Ann. Phys. 330(3), 377–445 (1908). [CrossRef]  

32. J. J. Faran Jr, “Sound scattering by solid cylinders and spheres,” J. Acoust. Soc. Am. 23(4), 405–418 (1951). [CrossRef]  

33. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University Press, 1999), 7th ed.

34. M. P. McLaughlin, “Compendium of common probability distributions,” (2016).

35. J. J. Carroll-Nellenback, R. J. White, R. W. Wood, and K. J. Parker, “Liver backscatter and the hepatic vasculature’s autocorrelation function,” Acoustics 2(1), 3–12 (2020). [CrossRef]  

36. D. S. Sivia and J. Skilling, Data Analysis: A Bayesian Tutorial (Oxford Science Publications, 2006), 2nd ed.

37. A. Clauset, C. R. Shalizi, and M. E. J. Newman, “Power-law distributions in empirical data,” SIAM Rev. 51(4), 661–703 (2009). [CrossRef]  

38. K. A. Vermeer, J. Mo, J. J. A. Weda, H. G. Lemij, and J. F. de Boer, “Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography,” Biomed. Opt. Express 5(1), 322–337 (2014). [CrossRef]  

39. Q. H. Vuong, “Likelihood ratio tests for model selection and non-nested hypotheses,” Econometrica 57(2), 307–333 (1989). [CrossRef]  

40. H. Akaike, “A new look at the statistical model identification,” IEEE Trans. Autom. Control 19(6), 716–723 (1974). [CrossRef]  

41. A. Kholodnykh, I. Petrova, M. Motamedi, and R. Esenaliev, “Accurate measurement of total attenuation coefficient of thin tissue with optical coherence tomography,” IEEE J. Sel. Top. Quantum Electron. 9(2), 210–221 (2003). [CrossRef]  

42. W. Yuan, C. Kut, W. Liang, and X. Li, “Robust and fast characterization of OCT-based optical attenuation using a novel frequency-domain algorithm for brain cancer detection,” Sci. Rep. 7, 44909 (2017). [CrossRef]  

43. S. Liu, Y. Sotomi, J. Eggermont, G. Nakazawa, S. Torii, T. Ijichi, Y. Onuma, P. W. Serruys, B. P. F. Lelieveldt, and J. Dijkstra, “Tissue characterization with depth-resolved attenuation coefficient and backscatter term in intravascular optical coherence tomography images,” J. Biomed. Opt. 22(9), 1–16 (2017). [CrossRef]  

44. S. Chang and A. K. Bowden, “Review of methods and applications of attenuation coefficient measurements with optical coherence tomography,” J. Biomed. Opt. 24(9), 090901 (2019). [CrossRef]  

45. P. Gong, M. Almasian, G. van Soest, D. M. de Bruin, T. G. van Leeuwen, D. D. Sampson, and D. J. Faber, “Parametric imaging of attenuation by optical coherence tomography: review of models, methods, and clinical translation,” J. Biomed. Opt. 25(4), 040901 (2020). [CrossRef]  

46. L. D. Fiske, M. C. G. Aalders, M. Almasian, T. G. van Leeuwen, A. K. Katsaggelos, O. Cossairt, and D. J. Faber, “Bayesian analysis of depth resolved OCT attenuation coefficients,” Sci. Rep. 11(1), 2263 (2021). [CrossRef]  

47. Y. Bromberg and H. Cao, “Generating non-rayleigh speckles with tailored intensity statistics,” Phys. Rev. Lett. 112(21), 213904 (2014). [CrossRef]  

48. L. Risser, F. Plouraboué, A. Steyer, P. Cloetens, G. Le Duc, and C. Fonta, “From homogeneous to fractal normal and tumorous microvascular networks in the brain,” J. Cereb. Blood Flow Metab. 27(2), 293–303 (2007). [CrossRef]  

49. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]  

50. I. Ben, Y. Y. Layosh, and E. Granot, “Study of a simple model for the transition between the ballistic and the diffusive regimes in diffusive media,” J. Biomed. Opt. 21(6), 066004 (2016). [CrossRef]  

Supplementary Material (1)

NameDescription
Supplement 1       Supplemental Info

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. Correspondence should be addressed to the corresponding author.

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 (5)

Fig. 1.
Fig. 1. Determining an ROI metric. (a) Simulated B-mode image using a Burr type XII distribution ($d=10,\,b=3$) with attenuation and shadowing. (b) Depth and lateral profiles based on spatial integration laterally and along depth, respectively. (c) Multiple ROIs and MLE fitting were used to calculate $\hat {b}$ as a function of maximum percent change in profiles, $\Delta$.
Fig. 2.
Fig. 2. Mouse brain with cranial window. (a) Single unfiltered and unaveraged B-mode frame with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles demonstrating that the maximum percent change $\Delta$ does not exceed 50%. (c) MLE fit to the Burr type XII distribution in a linear scale with linear binning of amplitude histogram data. (d) Same as (c) but in a log-log scale with logarithmic binning for visualization of tail behavior. (e) MLE fit to the Lomax distribution in a linear scale with linear binning of intensity histogram data. (f) Same as (d) but in a log-log scale with logarithmic binning to visualize tail behavior. (g) MLE fit to the generalized logistic distribution in a linear scale with linear binning of log amplitude histogram data. (h) Same as (g) but in log-log scale with logarithmic binning. All MLE fits result in an exponent parameter of $\hat {b}=6.06$ (95% Confidence Interval: $[5.83,\,6.28]$) for mouse brain.
Fig. 3.
Fig. 3. Excised mouse liver in phosphate-buffered saline. (a) Single unfiltered and unaveraged B-mode frame with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles demonstrating that the maximum percent change $\Delta$ does not exceed 50%. (c) MLE fit to the Burr type XII distribution in a linear scale with linear binning of amplitude histogram data. (d) Same as (c) but in a log-log scale with logarithmic binning for visualization of tail behavior. (e) MLE fit to the Lomax distribution in a linear scale with linear binning of intensity histogram data. (f) Same as (d) but in a log-log scale with logarithmic binning to visualize tail behavior. (g) MLE fit to the generalized logistic distribution in a linear scale with linear binning of log amplitude histogram data. (h) Same as (g) but in log-log scale with logarithmic binning. All MLE fits result in an exponent parameter of $\hat {b}=5.95$ (95% Confidence Interval: $[5.70,\,6.21]$) for mouse liver.
Fig. 4.
Fig. 4. Pig cornea. (a) Single unfiltered and unaveraged B-mode image with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles that do not exceed a $\Delta$ of 50%. (c) MLE fits to normalized amplitude with the Burr, Rayleigh, K, and Gamma distributions on a linear scale with linear binning. (d) Same as part (c) except on a log-log scale with logarithmic binning to visualize tail behavior. (e) MLE fits to normalized intensity with the Lomax, Exponential, K, and Gamma distributions on a linear scale with linear binning. (f) Same as part (e) except on a log-log scale with logarithmic binning to visualize tail behavior.
Fig. 5.
Fig. 5. Human hand (back side). (a) Single unfiltered and unaveraged B-mode image with the ROI shaded green. (b) The ROI’s normalized spatial integration profiles that do not exceed a $\Delta$ of 50%. (c) MLE fits to normalized amplitude with the Burr, Rayleigh, K, and Gamma distributions on a linear scale with linear binning. (d) Same as part (c) except on a log-log scale with logarithmic binning to visualize tail behavior. (e) MLE fits to normalized intensity with the Lomax, Exponential, K, and Gamma distributions on a linear scale with linear binning. (f) Same as part (e) except on a log-log scale with logarithmic binning to visualize tail behavior.

Tables (1)

Tables Icon

Table 1. Estimated exponent parameter values along with the 95% confidence interval for all samples. The average attenuation coefficient μ ^ avg in the ROI is also reported with the standard deviation. CI - confidence interval.

Equations (10)

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

p ( a ) = b 1 a min ( a a min ) b
I ( a ) = I 0 ( a a min )
p ( x ; d , b ) = 2 x ( b 1 ) d 2 [ ( x d ) 2 + 1 ] b , x > 0
p ( x ; d , b ) = ( b 1 ) d b 1 ( x + d ) b , x > 0
p ( y ; d , b ) = 2 ( b 1 ) exp ( 2 y ) d 2 [ exp ( 2 y ) d 2 + 1 ] b , < y <
p ( A ) = 2 A A 2 exp ( A 2 A 2 ) , A > 0 p ( I ) = 1 I exp ( I I ) , I > 0
p ( A ; α ) = 4 Γ ( α ) α A 2 ( α A 2 A 2 ) α / 2 K α 1 ( 2 α A 2 A 2 ) , A > 0 p ( I ; α ) = 2 Γ ( α ) α I I ( α I I ) α / 2 K α 1 ( 2 α I I ) , I > 0
p ( A ; α , β ) = 1 Γ ( α ) β α A α 1 exp ( β A ) , A > 0 p ( I ; α , β ) = 1 2 Γ ( α ) β α I α 2 1 exp ( β I ) , I > 0
μ ^ = 1 2 d z log ( 1 + I [ i ] i + 1 I [ i ] )
AIC = 2 k 2 log ( L ^ )
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.