Local Burr distribution estimator for speckle statistics

Speckle statistics in ultrasound and optical coherence tomography have been studied using various distributions, including the Rayleigh, the K, and the more recently proposed Burr distribution. In this paper, we expand on the utility of the Burr distribution by first validating its theoretical framework with numerical simulations and then introducing a new local estimator to characterize sample tissues of liver, brain, and skin using optical coherence tomography. The spatially local estimates of the Burr distribution’s power-law or exponent parameter enable a new type of parametric image. The simulation and experimental results confirm the potential for various applications of the Burr distribution in both basic science and clinical realms. © 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement


Introduction
The random interference pattern from coherent illumination known as speckle has been intensively studied for over a century. In the earliest studies, glass prisms were used to separate narrow color bands to investigate the phenomenon [1]. In the later 20 th century, optical speckle research was invigorated by the dissemination of radar and lasers [2][3][4][5][6][7][8][9]. In medical ultrasound, Rayleigh's early formulations were applied to speckle patterns that are seen in B-scan images [10]. Subsequent work applied a range of well-established statistical models for the backscattered echoes from tissues [11,12].
An important use of speckle metrics in medicine is the differentiation of normal versus diseased tissues. Newer models have been proposed [26,27] and employed to different clinical objectives [28][29][30][31]. A recent review of speckle in acoustics including the major mathematical formulations for probability density functions (PDFs) is presented in Stanton et al. [32]. They stress the influence of the scattering shapes and their distributions together with the beampattern effects.
Separately, in the field of optical coherence tomography (OCT) imaging of tissues, a parallel set of models have been employed both for speckle elimination and for physical modeling or biomarker characterization. Among these studies, OCT speckle analysis has generally emphasized utilization of the Rayleigh distribution [33][34][35][36][37], the K distribution [38], the gamma distribution [39,40], and comparisons between common models and a non-parametric approach [41].
More recently, a new speckle model of medical ultrasound based on power-law distributions originating in tissue vasculature was proposed [42][43][44][45]. The mathematics capture the behavior of a fractal branching set of cylinders modeled as a set of weak scatterers governed by a power-law distribution of scattering elements' radii. This framework breaks from earlier work with important distinctions. First, the classical assumption of many identical scatterers interacting so as to justify the central limit theorem is eliminated. The new model considers one or very few individual scatterers, and those are distributed in size with a power-law distribution. This led to a novel result -that ultrasound speckle distributions from vascularized tissue were effectively modeled by a two-parameter Burr distribution, originally formulated in the context of general functions for statistics research [46], which possesses a long power-law tail.
The Burr distribution enables a a new set of parameters and subsequent analyses that can be applied to characterize tissue imaged using OCT or ultrasound. This paper recaps the basic theory and then develops a relationship for spatially local estimates of the power-law parameter based on moments of the Burr distribution. Also, the effect of the number of scatterers included within a coherence volume is addressed by examining the complex addition of an increasing number of phasors. The Rayleigh distribution is used as a common comparison given its historical role in characterizing speckle in OCT and ultrasound [8,10,37]. Applications of these metrics to tissues are shown and assessed for their utility in tissue characterization.

Sufficient conditions for generating Burr speckle
Only two generalized conditions are sufficient to produce backscatter amplitudes that follow a Burr distribution. The first is a power-law distribution of weak scatterers, and the second is a linear increase of intensity with size or scale across some extended scale. Consistent with many fractal distributions [47] found in the natural world [48], we assume that the probability P of encountering a scatterer of dimension a within the ensemble is given by where a min is some minimum value of dimension, and b is the key power-law parameter governing the distribution of scattering size. In natural 3D fractals, b is typically found to be between 2 and 3 [47,49]. It is assumed that the scale of this distribution is valid from some minimum, small radii compared to the wavenumber k where k · a min ≪ 1 up to the point spread function of the imaging system. At larger sizes, targets begin to appear as discrete objects and are frequently excluded from the region of interest (ROI) analyzed as "speckle." Next, we submit that the backscatter intensity returning from a plane wave of unit amplitude encountering a scatterer of dimension a is assumed to be a linear function of a, above some minimum level set by the system and by Rayleigh scattering lower limits. This can simply be considered as a first order approximation of more complicated relationships that trend upward with increasing size, as shown by Parker [42]. Consequently, the amplitude A is a square root mapping function given as where A 0 is a constant scale factor. Given Eqs. (1) and (2), we can apply the probability transformation mapping rule [32,45,50] that directly yields the two-parameter Burr distribution [45]: The mean value of A is given as where Γ(·) is the gamma function. Thus, the mean or expected value of this distribution is proportional to the scale factor d; for the special case of b = 2.5, which is within the typical range for tissue, the mean is exactly equal to d. Furthermore, by direct integration, the ratio of the first moment (mean) squared and the non-centered second moment can be derived as Equation (5) is only dependent on b and increases monotonically towards a Rayleigh asymptote for b growing large. Figure 2 shows the monotonic increase in the estimator of Eq. (5) with increasing b from 2 to 8. The constant upper line in red is the same ratio of the first moment (mean) squared divided by the second moment for a Rayleigh distribution that can be shown to be equal to π/4 or approximately 0.8. Note this is not the same as the well-known Rayleigh signal-to-noise ratio of 1.91 because we are employing the non-centered second moment in the denominator instead of the traditional variance centered around the mean. Equation (5) in the limit of b → ∞ approaches the limit of π/4 exactly, thus suggesting an approach to the central limit theorem and the Rayleigh distribution. In practice, given a local region of speckle, we can easily calculate the first and second moments of the sampled data and solve Eq. (5) for b, which serves as our estimated power-law parameterb. Thus, the use of this ratio, Eq. (5), as a local estimator is proposed in this paper.

Additional factors: summation of several scatterers within the resolution cell
When multiple scatterers are located within the transmitted coherence volume, then a complex summation of their amplitudes is formed. When the phasors are characterized by PDFs, the summation of independent random variables is calculated using convolution or related techniques (see section IV of Stanton et al. [32]). It has been shown that the summation of Burr distributed random variables leads to series solutions [51][52][53][54][55][56][57]. However, for our case we need to add phasors The dependence of the ratio of the first moment (mean) squared divided by the non-centered second moment for the Burr distribution as a function of the Burr parameter b is shown in the blue curve. This ratio is a constant of π/4 for the Rayleigh distribution, shown in red. Note that as the power-law parameter b increases, the ratio asymptotically approaches the value found from a Rayleigh distribution. of the form |A| cos (θ) where A is a Burr distributed amplitude and θ is a uniformly distributed angle with 0 ≤ θ<2π. An explicit form of the expression can be found in [45], but a brief summary is given here. The PDF for |A| cos (θ) is derived from the product distribution rule for the multiplication of two random variables [32,50], which involves an integration. The following Pearson type VII distribution is thus obtained as where A ′ = |A| cos(θ). The Pearson type VII distribution is shown in Fig. 3. Next, the summation of two phasors is found by the convolution of two Pearson PDFs. Closed form analytic solutions can be found for specific half-integer values of b, so for a practical example [58] if b = 3.5 and d = 1, we derive the PDF for the sum of two Burr phasors as a polynomial ratio: With a summation of three phasors (i.e., an additional convolution of the PDFs), we obtain a higher order polynomial given as Following Rayleigh's derivations and arguments from August, 1880 [59], we can assume that the axes defining the phasor angles are arbitrary and consistent with an isotropic distribution of phasors expected over angle θ. As will be seen in the next sections, these PDFs will produce envelope distributions that strongly resemble Burr distributions over the peak area of the distribution; however, the estimated b will increase with an increasing number of convolutions. Nonetheless, at the extreme tail of rare, high amplitudes, we note that both Eqs. (7) and (8) asymptotically decay to ∼ A −6 , which corresponds to 2b − 1 = 6 in this case. In general, more phasor additions will yield PDFs with higher order polynomials commensurate with increasing estimates of Burr b as a curve fit result; however, the high amplitude (rare event) asymptotes remain limited.
Practically speaking, this implies that for lower resolution systems (with larger coherence volumes) the estimate of Burr parameter, denoted asb, will exceed the inherent b of Eq. (1) due to phasor addition. Historically, it has been widely assumed that as the number of scatterers becomes large, the result will be defined by the central limit theorem, leading to a Rayleigh PDF in amplitude [60]. In our case, beginning with the power law of Eq. (1), we anticipate only a few scatterers per sample volume within a high-resolution imaging system. Furthermore, the convergence of heavy-tailed power-law distributions is notoriously slow, so invoking the central limit theorem is not justified unless some special conditions are present, such as a broadly unfocused beam capturing a large volume of scatterers. These issues will be further illustrated in the Results section 4.

Influence of beampatterns
Equations (1) and (2) leading to the Burr PDF strictly apply to an ensemble of individual echoes or backscatter encountered by a uniform pulse or beam. In earlier work on small scatterers located randomly within an unfocused beampattern, the echo was modeled as a simple product of the random incident beam amplitude resulting from a random location within a prescribed beampattern multiplied by a random scattering transfer function [32]. The product rule for probability functions is used to derive the final overall distribution of echoes. However, when we now consider a multiscale distribution which included larger scatterers within a tightly focused beam, this assumption becomes questionable, and so other forms of simulations or models are required [43,61,62].

Numerical simulations
All numerical simulations were performed in MATLAB R2020b (MathWorks, Natick, Massachusetts, USA). Simulations were performed to compare Rayleigh and Burr distributions. The parameter m represents the number of phasors to be added, whereas n is the number of realizations. For each Rayleigh realization, m phasors of unit amplitude and random angles uniformly distributed between 0 and 2π were summed. The resultant summed vectors' amplitudes over n realizations form the histogram data. Upon normalization by the root mean square (RMS), the simulated data was then visually compared with the Rayleigh distribution.
For the Burr simulation, a power-law distribution with parameters b and a min was first created according to Eq. (1). Equation (2) was then calculated and the resulting distribution was applied as the phasors' amplitudes. The phasors' angles, θ, were still random and uniformly distributed between 0 and 2π. Upon summation of m phasors for all n realizations, the resultant summed amplitudes were normalized by dividing by the RMS, and subsequently fitted to the two-parameter Burr distribution as prescribed by Eq. (3) using the maximum likelihood estimation (MLE) method. The built-in MLE function in MATLAB was used and Eq. (3) was the input for the custom distribution.
Simulations were performed by varying the parameter m = [1,2,4,8,16], while setting b = 3, a min = 1, and n = 1, 000, 000. The choice for b is arbitrary, but we chose the numerical value such that it is within our recently tabulated estimated values from tissues [58]. We chose a million sample points, which based on our previous analyses, is a sufficient ensemble to capture accurate behavior without excessively increasing computational time. In addition, the estimated power-law parameterb is measured as a function of varying m phasors for b = [2, 2.5, 3] using both the MLE method and based on Eq. (5). These are respectively denoted asb MLE andb EE (where EE denotes equation estimate). The estimated parameterd is estimated only using MLE.. Finally, the same simulations were used to verify the Pearson distribution, or Eq. (6), from Section 2.2 by incorporating the factor of cos(θ) in the analysis.

Experimental setup
A swept source optical coherence tomography (SS-OCT) system was used to scan samples of calf liver, mouse brain, and human skin. It was implemented with a swept source laser (HSL-2100-WR, Santec, Aichi, Japan) with a center wavelength of 1310 nm and bandwidth of 170 nm. The lateral resolution is 20 µm and the axial resolution is 8 µm in air. The maximum sensitivity of the system was measured to be approximately 120 dB. The imaging depth was measured to be 5 mm in air. The SS-OCT system was controlled with LabVIEW (Version 14, National Instruments, Austin, Texas, USA). The scans (single frame, unaveraged, unfiltered, uncorrected) were performed with 500 A-lines and the focal plane was set to approximately 0.2 mm below the sample surface. The mouse brain was scanned in vivo using a cranial window that enables OCT imaging. The mouse study was performed in accordance with experimental protocols approved by the University of Rochester Committee on Animal Resources.

Burr local estimation
For each scan, a global estimate in the ROI is obtained using MLE by implementing the procedure outlined in Ge, et al. [58]. This procedure includes a test for uniformity of amplitudes within the ROI so as to avoid shadows or dissimilar tissues. Another global estimate using Eq. (5) is obtained for direct comparison. Finally, a local map of the Burr power-law parameter is obtained by applying the estimator of Eq. (5) within moving a rectangular window of 100 × 80 pixels across the ROI with intervals of 5 × 4 pixels. While the selection of this window size and interval is arbitrary, we find the windows are adequately sized so that local estimates are accurate, while resolution is also preserved. Linear interpolation is performed after the estimates are obtained to fill out the ROI, and no other filters are applied.

Numerical simulations
First, we examine the background issue of the amplitude of the sum of random phasors, comparing the classical case of uniform amplitude phasors with Burr distributed phasors. In the series shown in Fig. 4, it is apparent that as the number of phasors increases, the Rayleigh distribution is approximated with eight identical phasors or more. The Burr estimators accurately estimate b with a single phasor, and overestimates b as the number of phasors increases. We report that b EE is always within the 95% confidence interval obtained withb MLE for each case. The scale parameterd also increases as more phasors are summed. With a large number of phasors, the Burr estimate resembles that of the Rayleigh distribution. The last subfigure of panel 4(b) is shown in Fig. 5, on a log-log scale to demonstrate the differences among the data points, the Burr PDF, and the Rayleigh PDF. In Fig. 6(a), it is appreciated that at m = 1,b matches b. As the number of phasors increases, the estimatedb increasingly overestimates b. A similar trend is noted withd in Fig. 6(b). Figure 7 demonstrates the validity of the Pearson distribution (i.e., Eq. (6)), as explained in Section 2.2. With the phasor magnitudes incorporating the cos(θ) random angle, the Pearson distribution can also be used to estimate the Burr parameter b.

Burr local estimation in tissue
In the following examples, we compare results from global estimates ofb obtained from the entire ROI, against the local estimates within the sliding window used to makeb images. The global MLE parameter for the ROI of the OCT liver scan shown in Fig. 8(a) isb MLE = 5.24.   The global estimate using Eq. (5) isb EE = 5.20. The local map using Eq. (5) on local, sliding windows is shown in Fig. 8(b), which shows strong homogeneity in theb local estimate. The global MLE parameter for the ROI of the OCT mouse brain scan shown in Fig. 9(a) isb MLE = 4.93. The global estimate using Eq. (5) isb EE = 4.94. The local map is shown in Fig. 9(b), which shows homogeneity, similar to the case with liver. The global MLE parameter for the ROI of the OCT skin scan shown in Fig. 10(a) isb MLE = 3.12. The global estimate using Eq. (5) isb EE = 3.14. The local map is shown in Fig. 10(b), which shows a visual correlation between the layers of the skin and theb local estimate.

Discussion and conclusion
The Burr distribution emerges from two sufficient conditions. First, the probability that a pulse or beam encounters a scatterer of some dimension is assumed to be determined by a power law. Secondly, the backscatter intensity is assumed as a first order approximation to increase linearly with size within the sub-resolved scale. By applying the basic rule of probability transformation, these two assumptions generate speckle PDFs formed by a two-parameter Burr distribution. The particular Burr distribution we use is governed by a scale parameter d and a power-law exponent b. This Burr model of speckle amplitudes is appropriate for scatterer distributions where multi-scale or fractal distribution of structures is present and is now found to be a useful framework in both medical ultrasound and OCT [41,58].
Numerical simulations using phasor addition are utilized to demonstrate how the two-parameter Burr PDF compares with the classical Rayleigh PDF. As the number of phasors increases, the Rayleigh simulations approach the Rayleigh distribution, which attests to its underlying assumption of a large number of scatterers. For the Burr simulations, with a single phasor, accurate powerlaw parameters can be estimated. As the number of phasors increases, the parameters are overestimated, and the simulation distribution eventually converges to the Rayleigh PDF.
A new local estimator using the basis of Eq. (5) is also proposed and demonstrated with three sample OCT scans. As seen in the liver and brain scans, the estimated map ofb is homogeneous. In the skin scan, the estimated map ofb is heterogeneous, and closely follows the layers of skin, also distinctly outlined in the B-mode image. These three scans demonstrate the potential utility of a local Burr estimator to characterize tissue. However, further studies are needed to investigate its performance and limitations with respect to different ranges of the mean index of refraction, multiple scattering, beampattern effects, and window optimization. We have not emphasized the parameterd which is related to the average amplitude of the speckle and is strongly dependent on the gain settings of the imaging system. It is possible to construct methods to normalize the system parameters and extract a tissue-based estimate of d related to optical scattering. These studies are left to future research.
In summary, the utility of the Burr distribution has been expanded upon with numerical simulations to validate theory as well as experiments to demonstrate its potential in identifying and characterizing local tissue structures. New and alternative estimators proposed can be further refined with various image processing techniques, which would allow for translation to direct biological and clinical applications.