An Improved Parameter Estimator of the Homodyned K Distribution Based on the Maximum Likelihood Method for Ultrasound Tissue Characterization

Abstract

The homodyned K distribution (HK) can generally describe the ultrasound backscatter envelope statistics distribution with parameters that have specific physical meaning. However, creating robust and reliable HK parameter estimates remains a crucial concern. The maximum likelihood estimator (MLE) usually yields a small variance and bias in parameter estimation. Thus, two recent studies have attempted to use MLE for parameter estimation of HK distribution. However, some of the statements in these studies are not fully justified and they may hinder the application of parameter estimation of HK distribution based on MLE. In this study, we propose a new parameter estimator for the HK distribution based on the MLE (i.e., MLE1), which overcomes the disadvantages of conventional MLE of HK distribution. The MLE1 was compared with other estimators, such as XU estimator (an estimation method based on the first moment of the intensity and tow log-moments) and ANN estimator (an estimation method based on artificial neural networks). We showed that the estimations of parameters α and k are the best overall (in terms of the relative bias, normalized standard deviation, and relative root mean squared errors) using the proposed MLE1 compared with the others based on the simulated data when the sample size was N = 1000. Moreover, we assessed the usefulness of the proposed MLE1 when the number of scatterers per resolution cell was high (i.e., α up to 80) and when the sample size was small (i.e., N = 100), and we found a satisfactory result. Tests on simulated ultrasound images based on Field II were performed and the results confirmed that the proposed MLE1 is feasible and reliable for the parameter estimation from the ultrasonic envelope signal. Therefore, the proposed MLE1 can accurately estimate the HK parameters with lower uncertainty, which presents a potential practical value for further ultrasonic applications.

Keywords

homodyned K distribution statistical parameter estimation maximum likelihood estimation quantitative ultrasound ultrasonic tissue characterization

Introduction

Biopsies often serve as the gold standard for diagnosing certain diseases. However, a biopsy is associated with invasive procedures, sampling errors, and time and cost constraints. It is unnecessary particularly in the case of misdiagnosis or ambiguity.^1
-3 Medical imaging, particularly ultrasound imaging, has emerged as a noninvasive alternative to biopsies. Compared to other medical imaging modalities, ultrasound imaging incurs a lower cost, is user-friendly, and uses non-ionizing radiation. Therefore, it has gained widespread acceptance and has been widely used in routine diagnostics.^2,4,5 Ultrasound echo signals are generated by a process in which a transducer transmits ultrasound waves into a tissue being examined and receives the backscattered signal, which is caused by the interaction between the tissue and the incident waves.⁴ Conventional B-mode ultrasound images are constructed from beam-formed radiofrequency (RF) echo signals after they have been subjected to scan conversion (resampling onto a rectangular grid), logarithmic compression, and other proprietary nonlinear mappings built into the ultrasound machine.⁶ These processes are necessary to make raw ultrasound data readable to the human eye. However, they inevitably lose some useful properties that are theoretically derived from the RF signal. In addition, conventional B-mode images lack functional information and suffer from the drawbacks of being qualitative and operator-dependent.^2,5
-10 Therefore, quantitative ultrasound (QUS) techniques, a complement to conventional B-mode ultrasound, have been proposed. In general, QUS refers to any method that extracts quantitative tissue information from ultrasound backscattered signals for ultrasound tissue characterization. Moreover, QUS techniques consist of spectral-based parameterization of ultrasound signals, flow estimation through Doppler, tissue elastography techniques, shear wave imaging, and envelope statistics.^4,11
-18

Ultrasound backscatter envelope statistics, both statistical model-based and non-model-based, is an important group of QUS techniques.² It extracts statistical information related to tissue microstructures from the envelope of backscattered RF signals and has been applied to characterize different types of tissue. In particular, statistical model-based ultrasound envelope statistics parametric imaging techniques, in contrast to non-model-based techniques, such as ultrasound kurtosis imaging and ultrasound entropy imaging, have received broad interest in the past three decades.^2,19
-23 Various statistical models have been proposed to estimate the information regarding sub-resolution material properties.^4,24 Among them, the Rayleigh distribution is the most basic and provides the least amount of information. The Rayleigh distribution arises when numerous nearly identical and randomly located scatterers contribute to the echo signal.¹¹ The K distribution is another model that is considerably more general than the Rayleigh distribution. It can provide information on the scatter density, variation in the scattering amplitudes within the range cell, and mean scattering amplitude from the two parameters of the K distribution, namely, the effective number of scatterers and the scaling parameter.^25
-27 The K distribution approaches the Rayleigh distribution under conditions of high values of the effective number of scatterers.^4,11,24 The K distribution, although better than the Rayleigh distribution, is not general enough to account for the existence of periodic scatterers. This flaw can be addressed by adopting the Nakagami distribution or homodyned K (HK) distribution.^25,26,28 Nakagami distribution, a good approximation of the HK distribution, has been the most frequently adopted model in tissue characterization because of its simplicity, even though the HK distribution is more versatile.^1,4,26 In recent years, however, comparative studies have reported that the HK parameter, α, is more appropriate than the Nakagami parameter, m, to describe media with a high number of scatterers per resolution cell.^4,7 In addition, the parameters of the HK distribution contain more information regarding the tissue microstructure than those of the Nakagami distribution.^1,4 Therefore, with the development of effective parameter estimators, the HK model has been used in ultrasound parametric imaging for various medical applications, such as classification of breast lesions, temperature monitoring during focused ultrasound treatment, and hepatic steatosis evaluation.^{2,7,21,29
-31}

To estimate the parameters of the statistical distribution models, methods based on maximum likelihood estimator (MLE), moments, or artificial neural networks (ANNs) are commonly adopted. For the HK distribution, Dutt and Greenleaf³² first proposed a method based on the moments of the intensity (i.e., MI estimator). A method based on the first three moments of the amplitude was used to estimate the HK distribution in 1995 (i.e., MA estimator), but a closed-form expression of the amplitude moments could not be found until an explicit expression of an arbitrary moment of the amplitude was given by Hruska and Oelze.³⁴ Moreover, improved estimators for the HK distribution, which are based on the integer or fractional moments of envelope signals and ANNs, have been proposed. They include (i) the estimator based on the signal-to-noise (SNR) of fractional order moments of the amplitude³³; (ii) the level-curve estimator based on the SNR, skewness, and kurtosis of the fractional order moments of the amplitude, namely the RSK estimator³⁴; (iii) the estimator based on the mean intensity and two log-moments, namely the XU estimator³⁵; and (iv) the estimator based on an ANN, namely the ANN estimator.³⁶

MLE usually has desirable asymptotic properties.^37,38 Additionally, comparative studies have presented parameter estimators based on MLE or maximum a posteriori estimation that have a small normalized mean squared error in the case of the Rice, K, or Nakagami distribution.⁴ Thus, as a classical parameter estimation method, the MLE has been used to estimate the parameters of the HK distribution. Hu et al.³⁹ first used MLE to estimate parameters of HK distribution to assess the physical meanings of parameters and goodness of fit of HK distribution for ultrasound envelope data from scatterer distributions with wide-varying spatial organizations. They used a high-order global adaptive Gauss-Kronrod quadrature method to calculate the probability density function (PDF) of the HK distribution, and then used the Newton-Raphson method to numerically solve the log-likelihood equation. Moreover, Wang et al.⁴⁰ proposed a similar MLE for parameter estimation of HK distribution but adapted a cloud adaptive particle swarm optimization (CAPSO) to calculate the PDF of the HK distribution. However, these studies have significant limitations, and some of the statements are misleading. We will further discuss them in detail in Section 2.5.

MLE can be used to develop an accurate and stable parameter estimation method for the HK distribution, however conventional MLE is flawed. Therefore, there is an urgent need to propose an improved MLE which can overcome the shortcomings of conventional MLE and facilitate the application of parameter estimator for HK distribution based on MLE in ultrasound tissue characterization. In addition, the sample size is an important factor affecting the performance of the HK distribution parameter estimator, so we explored the dependence of the proposed MLE1 performance on the sample size. Furthermore, inspired by the work in Cristea et al.,⁷ we further explored the usefulness of the proposed MLE1 for characterization of both sparse and concentrated media. The work of this study is two folds. First, we proposed a novel parameter estimation method of the HK distribution based on the MLE (MLE1) after analyzing the deficiencies of conventional MLE of HK distribution. Second, we evaluated the feasibility and performance of the proposed MLE1 through simulation experiments based on independent and identically distributed (i.i.d.) samples from the HK distribution and ultrasound images.

The remainder of this paper is organized as follows. In Section 2, the theoretical background of parameter estimation of the HK distribution is presented. In Section 3, the proposed estimation method is explained, and computer simulations for sample sets of the HK distribution and ultrasound images are described. The simulation results and discussions are presented in Section 4, and the conclusions are provided in Section 5.

Theoretical Background

Homodyned K Distribution

The PDF of the HK distribution is significantly complicated and does not have a closed form. In 1987, the expression of HK distribution was introduced by Jakeman and Tough²⁷ in terms of an improper integral. It (two-dimensional) is defined as follows:

P_{H K} (A | ε, σ^{2}, α) = A \int_{0}^{\infty} u J_{0} (u ε) J_{0} (u A) {(1 + \frac{u^{2} σ^{2}}{2})}^{- α} d u,

(1)

where A is the amplitude of the backscattered ultrasound envelope, $J_{0} (\cdot)$ denotes the Bessel function of the first kind of order 0, and α, ε, and σ are the three parameters of the HK distribution. In particular, α is the scatterer clustering parameter that is related to the density of the scatterers and the homogeneity of the scatterers. The mean intensity of the backscattered ultrasound can be expressed as

\bar{I} = ε^{2} + 2 α σ^{2},

(2)

where $ε^{2}$ and $2 α σ^{2}$ are the coherent and diffuse signal components, respectively. Additionally, the coherent-to-diffuse signal ratio is equal to

k = \frac{ε}{σ \sqrt{α}} .

(3)

The PDF of HK distribution can also be defined in a compound form.²⁷

P_{H K} (A | ε, σ^{2}, α) = \int_{0}^{\infty} R_{R i} (A | ε, σ^{2} w) ς (w | α, 1) d w,

(4)

R_{R i} (A | ε, σ_{R}^{2}) = \frac{A}{σ_{R}^{2}} I_{0} (\frac{ε}{σ_{R}^{2}} A) \exp (- \frac{(ε^{2} + A^{2})}{2 σ_{R}^{2}}),

(5)

ς (w | α, 1) = w^{α - 1} \exp (- w) / Γ (α),

(6)

where $P_{R i}$ denotes the Rice distribution, $ς (w | α, 1)$ is the gamma distribution with mean and variance α, and $Γ (\cdot)$ denotes the Euler gamma function.²⁴

RSK Estimator

The RSK estimator estimates the parameters α and k of the HK distribution based on the SNR, skewness S, and kurtosis K from the backscattered ultrasound envelope. It can be defined as^33,34

\begin{array}{l} R_{v} = \frac{E [A^{v}]}{{(E [A^{2 v}] - E^{2} [A^{v}])}^{1 / 2}}; \\ S_{v} = \frac{E [A^{3 v}] - 3 E [A^{v}] E [A^{2 v}] + 2 E^{3} [A^{v}]}{{(E [A^{2 v}] - E^{2} [A^{v}])}^{3 / 2}}; and \\ K_{v} = \frac{\begin{array}{l} E [A^{4 v}] - 4 E [A^{v}] E [A^{3 v}] + 6 E [A^{2 v}] E^{2} [A^{v}] \\ - 3 E^{4} [A^{v}] \end{array}}{{(E [A^{2 v}] - E^{2} [A^{v}])}^{2}}, \end{array}

(7)

where A is the amplitude of the backscattered ultrasound envelope, v is the power of the amplitude, and $E [\cdot]$ is the expectation operator. The estimates of the SNR, skewness, and kurtosis for the backscattered ultrasound envelope signal can be used to obtain parameter estimates by finding the point in the (α, k) space where the L2-norm of the distances from the level curves is minimized.³⁴

XU Estimator

The estimation method of the HK distribution based on the mean intensity and X- and U-statistics can be presented as follows³⁵:

E [I] = \bar{I}; U_{H K} = U; X_{H K} = X,

where X and U are defined as

\begin{array}{l} X = E [I l o g I] / E [I] - E [l o g I]; and \\ U = E [\log I] - logE [I] . \end{array}

(9)

The expressions of $U_{H K}$ and $X_{H K}$ have been previously defined³⁵ in Destrempes et al. (i.e., equations (3.3) and (3.4)). With respect to numerical considerations, in practice, the three parameters of the HK distribution are estimated using the following algorithm:

\begin{array}{l} (ε, σ^{2}, α) = \arg \min | U_{H K} - U | \\ subject to μ = \bar{I} and X_{H K} = X and, α \leq α_{m a x} \end{array}

(10)

where $α_{m a x}$ is the upper bound of the cluster parameter α, and $\bar{I}$ is the mean intensity of the backscattered ultrasound envelope signal.

ANN Estimator

Recently, a parameter estimator of the HK distribution based on an ANN was proposed, which has a much faster computing speed and a smaller normalized standard deviation for parameters α and k compared with the RSK and XU estimators.³⁶ The architecture of the ANN is shown in Figure 1, which is a four-layer feed-forward backpropagation (BP) network. It consists of an input layer, two hidden layers, and an output layer. By using equations (7) and (9), SNR R, skewness S, kurtosis K, and X- and U-statistics can be extracted from the i.i.d. envelope signal samples of the HK distribution with respect to different combinations of α and k parameters, and can be input into the corresponding neurons in the input layer. The numbers of neurons in the two hidden layers were 30 and 10, respectively. The output layer is made up of two neurons, which correspond to the $\log_{10} α$ and k estimates. This estimator can be easily implemented using the MATLAB platform. More information can be found in Zhou et al.³⁶ in detail.

Conventional MLE for Parameter Estimation of HK Distribution

The conventional MLE for parameter estimation of HK distribution mainly consists of two parts. A numerical integration method is used to calculate the PDF of the HK distribution, and then a numerical calculation method is used to find the parameters that maximize the log-likelihood equation.

First, the authors in Hu et al.³⁹ adapted a numerical integration function quadgk in the MATLAB platform to calculate the PDF of the HK distribution. Whereas the authors in Wang et al.⁴⁰ considered that significant oscillation is produced by the Bessel function contained in the PDF of HK distribution, therefore they used an advanced CAPSO algorithm to calculate the PDF of the HK distribution. However, such a statement is not fully justified. As mentioned above, the PDF of the HK distribution has two expressions. In fact, for different PDF expressions of HK distribution, the properties of their integrands are different. It can be seen from Figure 2 that for equation (1), its integrand oscillates and does not converge when α is small. In this case, the numerical integration method cannot be used to calculate the integral. Moreover, when the value of α is large, its integrand converges quickly and the function curve is very steep. Therefore, if equation (1) is used to calculate the PDF of the HK distribution, the performance of the numerical integration method will have a great influence on the calculation result. On the contrary, the integrand of equation (4) is convergent, and the function curve is flat. We can conclude that it is more reasonable to use equation (4) instead of equation (1) to calculate the PDF of the HK distribution. Thus, the oscillation may appear in Hu et al.³⁹ which uses equation (1) rather than in Wang et al.⁴⁰ which uses equation (4). Moreover, the performance of the numerical integration method has slight effect on the calculation results when adapting equation (4) and we will further verify it through simulation experiments.

Second, both the article Hu et al.³⁹ and Wang et al.⁴⁰ use the Newton-Raphson method to find the critical point that maximizes the log-likelihood equation. And the authors believe that sole solution can be obtained based on the Newton-Raphson algorithm to solve the log-likelihood equation numerically.⁴⁰ However, this statement is true unless the log-likelihood equation has only one critical point in the parameter domain because Newton-Raphson method can only find one root near the initial value. Recall that the Maximum likelihood estimation can be defined as a point $\bar{θ}$ _* in the parameter domain at which (i) the gradient of the log-likelihood function $\log (L (\bar{θ}$ )) vanishes (i.e., $\bar{θ}$ _* is a critical point of the log-likelihood function); (ii) $\log (L (\bar{θ}$ _*)) is maximal among all critical points of the log-likelihood function (at least in some bounded domain).⁴ In fact, the equation $\nabla \log (L (\bar{θ}$ )) = 0 may have several solutions. Therefore, a locally optimal solution may be found by using the Newton-Raphson method (i.e., $\bar{θ}$ _i is a critical point, but it cannot maximize the log-likelihood equation). Worse still, according to the statement in article Mamou and Oelze,⁴ the equation $\nabla \log (L (\bar{θ}$ )) = 0 may have no solution. In this case, the Newton-Raphson method will not converge. Therefore, we need to further explore a method that can find a point that maximizes the log-likelihood equation (at least in some bounded domain) and always converges. In addition, relevant contents such as the range of the various functions and the setting of the initial values in an implementation have not been sufficiently discussed in the article Hu et al.³⁹ and Wang et al.⁴⁰

Generation of Independent and Identically Distributed Samples From Homodyned K Distributions

To evaluate the estimators, Monte Carlo simulations are usually performed with i.i.d. samples from the HK distribution, which are generated with known parameters. The i.i.d. samples can be obtained as^33,34

A_{i} = \sqrt{{(ε + X_{i} σ \sqrt{z_{i}})}^{2} + {(Y_{i} σ \sqrt{z_{i}})}^{2}},

(11)

where $A_{i}$ denotes an i.i.d. HK sample, $X_{i}$ and $Y_{i}$ are i.i.d. samples of the standard Gaussian distribution, and $z_{i}$ is a sample of the gamma distribution with shape parameter α and scale parameter unity. It is worth noting that if we multiply all the samples by the same constant, it does not change the size of k, because σ and ɛ are multiplied by the same constant and α does not change. This is a very important property, and we will refer to it as property 1 later.

Moreover, if we restrict the average intensity to 1, we can obtain the following expression by using equations (2), (3), and (11).³⁶

A_{i} = \sqrt{{(\frac{k}{\sqrt{k^{2} + 2}} + X_{i} \frac{\sqrt{z_{i} / α}}{\sqrt{k^{2} + 2}})}^{2} + {(Y_{i} \frac{\sqrt{z_{i} / α}}{\sqrt{k^{2} + 2}})}^{2}}

(12)

Materials and Methods

In this section, first, we proposed a novel parameter estimator for HK distribution based on the MLE, using a numerical integration method to approximate the PDF of the HK distribution and a non-gradient global optimization algorithm to find a critical point that maximizes the log-likelihood function in the given parameter domain. Second, we carried out a series of simulation experiments based on the i.i.d. samples of HK distribution with known parameters. Specifically, (1) we adopted the method in Destrempes et al.³⁵ to test the proposed MLE1 through Monte Carlo computer simulation in the domain $k \in {0.1, 0.2, \dots, 0.9, 1.0} and α \in {1, 2, \dots, 9, 10}$ , and compared the results with those obtained from MI, MA, RSK and XU estimators in Destrempes et al.³⁵ to verify the feasibility and accuracy of the proposed MLE1. At the same time, different numerical integration methods were used in the proposed MLE1 to explore the influence of numerical integration performance on parameter estimation. (2) Inspired by the work in Cristea et al.,⁷ we further evaluated the usefulness of the proposed MLE1 to characterize the media with a wide range of numbers of scatterers per resolution cell through computer simulation, and compared the simulation results with those of the XU estimator in Cristea et al.⁷ (3) Considering the sample size may have an important impact on the results of QUS techniques based on envelope statistics. Therefore, we conducted simulation experiments to explore the dependence of the performance of the proposed MLE1 on the sample size, and conducted a similar study on the recently proposed ANN estimator for comparison. Third, tests on two types of simulated ultrasound images based on Field II (one possesses only a diffuse component and the other possesses both diffuse and coherent components) were performed to evaluate the feasibility and reliability of proposed MLE1 for the parameter estimation of ultrasonic envelope signal. In addition, the sensitivity and robustness of using different estimators (i.e., RSK, XU, ANN, and MLE1) in ultrasonic tissue characterization were assessed by comparing the estimated parameter values and the corresponding differences with theoretical values under different number densities of scatterers per resolution cell.

Proposed MLE

The MLE is a generalized statistical method to estimate parameters, and it is defined as a critical point that maximizes the log-likelihood function³⁸:

L (\bar{θ}) = \sum_{i = 1}^{N} l o g P (A_{i} | \bar{θ}),

(13)

where $\bar{θ}$ represents the parameters to be estimated, ${A_{1}, A_{2}, \dots \dots ., A_{N}}$ is the sample data of size N, and $P (A_{i} | \bar{θ})$ is the PDF of the statistical distribution. Usually, to obtain the critical point, we need to compute the partial derivatives of the log-likelihood function and solve a system of equations. However, calculating the partial derivatives of the log-likelihood function are complex and time-consuming. Therefore, a pattern search optimization algorithm was used in this study to determine the critical point of the global maximum of log-likelihood function in a given parameter domain. Pattern search algorithm is one of the most popular direct search methods that are designed to solve derivative-free problems. Unlike optimization methods that require information about the gradient or higher derivatives of the objective function, a direct search algorithm searches a set of points around the current point, looking for one where the value of the objective function is lower than the current point.⁴¹ It can be used to consider the MLE problem as a bound, constrained nonlinear program, with the main objective of minimizing $- L (\bar{θ}) .$ The pattern search algorithm is a global optimization algorithm and can be easily implemented using the patternsearch function in the MATLAB global optimization toolbox.⁴¹ It is worth mentioning that the parameter settings used in the patternsearch function were the default values, except as stated otherwise.

The flow charts of the proposed MLE1 for estimating the parameters of the HK distribution are shown in Figures 3 to 5. First, we obtained the i.i.d. HK distribution samples A_i to be fitted, and initialized the values of the parameters to be estimated for their upper and lower boundaries, respectively. In addition, the termination condition (i.e., MaxIterations and MeshTolerance) for this function needed to be specified, and the MeshTolerance was set at 0.001 in this study. Next, we used the pattern search algorithm to repeatedly update and iterate over the estimated parameters to minimize $- L (\bar{θ})$ globally. To calculate the value of $- L (\bar{θ})$ , we adopted numerical integration methods to calculate the PDF of the HK distribution. In Figure 4, we used the traditional trapezoidal numerical integration method, where the upper limit of integration was set to 200, the lower limit was set to 0, and the number of integration nodes was set to 1000. In Figure 5, we use the adaptive numerical method based on CAPSO.^40,42,43 The step size of the numerical integration can adaptively change by minimizing the fitness function $f t = \sum_{j = 1}^{n} \frac{1}{2} (x_{j + 1} - x_{j}) (f (x_{j + 1}) - f (x_{j}))$ to obtain a more accurate result. The maximum number of iterations MaxIterations was set to 20. The number of integration nodes was set to 1000. The upper and lower limits of integration were the same as those in Figure 4. One can refer to Zhu⁴² and Qu and He⁴³ for more details on numerical integration based on CAPSO.

The following aspects are worth discussing with regard to the proposed MLE1 method. First, one may consider how to select the initial values of the parameters to be estimated and their upper and lower boundaries. In general, based on values estimated from ultrasound data, the range of parameter α is between 1 and 80, and between 0 and 1 for parameter k.⁷ Additionally, according to property 1 discussed in the previous section, we can limit the value range of the HK distribution samples. For example, if we multiply all the samples by $10 / \max (A_{i})$ , all the values of the samples are restricted to be between 0 and 10. The mean intensity $\bar{I}$ calculated from the HK distribution samples is less than or equal to 100. Therefore, we can estimate the range of parameters ɛ and σ based on equation (2). In the aforementioned case, the range of 0 to 10 is sufficient for parameters ɛ and σ. In practice, we can first select a large range of parameter domain for parameter estimation, and then reasonably narrow the parameter space according to the obtained estimated value to improve the efficiency and accuracy of the algorithm. As for the setting of the initial values of the parameters to be estimated, if the approximated ranges of the theoretical values are unavailable, we can set them as the upper boundary value of each parameter to be estimated. This method is very effective in practical applications.

Second, the upper limit of numerical integration must be chosen. According to the PDF of the HK distribution, it is known that the upper limit is supposed to be set to positive infinity in theory, but this is unrealistic in practice. Meanwhile, a larger upper limit corresponds to a large computational cost, which is time-consuming. Fortunately, this challenge can be dealt with well by using property 1. Based on the aforementioned example, when we limit the maximum value of all sample values to 10, the range of parameters ɛ and σ can be 0 to 10, and the range for parameter α is 0 to 80. In this parameter domain, we found that setting the upper limit of integration to 200 is sufficient to calculate the PDF well by comparing the PDF curve calculated by the numerical integration method with the histogram envelope of samples generated using equation (11) under the same parameters.

Finally, as the PDF of the HK distribution contains a Bessel function, it may lead to a very large value. Thus, there is a problem with the range of parameters to be estimated with respect to numerical considerations. Considering that the maximum positive floating point number of the MATLAB platform is $1.7977 \times 10^{308}$ , when the variable $\frac{ε A}{σ^{2} w}$ in the zero-order first-kind Bessel function is greater than 700, it causes the stored data to overflow. Therefore, we hope that the sample value, A_i, is not too large, and the values of the first few numerical integration nodes (i.e., w₁, w₂,. . .) are not too small. In practice, we can use property 1 to limit the size of the sample value and ignore the first few numerical integration nodes after weighing the feasibility and accuracy of numerical integration. Moreover, since pattern search algorithm can search for global solutions for problems that contain multiple maxima or minima, the proposed MLE1 can find a point of the global maximum of the log-likelihood in a given parameter domain.⁴¹

Simulations of Independent and Identically Distributed Samples From Homodyned K Distributions

To evaluate the performance of the proposed MLE1, we tested it through several Monte Carlo computer simulations by generating sets of i.i.d. samples of the HK distribution with known parameters. We adopted the same approach as in Destrempes et al.³⁵ to compare the reported experimental results. Sets of samples were generated in the domains $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ . In the subsequent simulations of i.i.d. samples from the HK distribution, the value of $σ^{2}$ was taken as 1, so that the value of $ε$ was set to $k \sqrt{α}$ . In the pattern search algorithm, the initial values of parameters ɛ $and α,$ which are to be estimated, were set to 4 and 10, respectively. Their lower boundaries were set to 0.01 and 0, and their upper boundaries were set to 4 and 10, respectively. For each combination of parameters, 1000 samples of size N = 1000 were simulated using equation (11). For the numerical integration method, the traditional trapezoidal numerical integration method and the numerical integration method based on CAPSO were used to explore the influence of the calculation accuracy of the PDF for the HK distribution on the parameter estimation. The results of parameter estimates using the MLE1 were compared with the experimental results using the MI, MA, XU, and RSK estimators in Destrempes et al.³⁵ For each estimation method, the relative biases (RB), $(E [\hat{k}] - k) / k$ and $(E [\hat{α}] - α) / α$ , the normalized standard deviations (NSD), $\sqrt{V a r [\hat{k}]} / k$ and $\sqrt{V a r [\hat{α}]} / α$ , and the relative root mean squared errors (RRMSE), defined by $\sqrt{E [{(\hat{k} - k)}^{2}}] / k$ and $\sqrt{E [{(\hat{α} - α)}^{2}}] / α$ , were calculated based on 1000 estimates of $\hat{k}$ and $\hat{α}$ with those parameters.

Previous simulation studies mostly focused on low numbers of scatterers per resolution cell, while it is meaningful to evaluate scattering properties from concentrated media. Inspired by the work of Cristea et al.,⁷ we assessed the usefulness of the proposed MLE1 for estimating concentrated media using i.i.d. samples from the HK distribution. Sets of samples were generated in the domains $k \in {0.0, 0.1, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 79, 80}$ . In the pattern search algorithm, the initial values of the parameters ɛ $and$ α, which are estimated, were set to 10 and 80, respectively. Their lower boundaries were set to 0.01 and 0, and their upper boundaries were set to 10 and 80, respectively. For each combination of parameters, 100 samples of size N = 1000 were simulated using equation (11). Here, the numerical integration method based on CAPSO was used to calculate the PDF of the HK distribution. To compare with the simulation results in Cristea et al.,⁷ we used the same evaluation metrics (i.e., standard deviations (SD) $\sqrt{V a r [\hat{k}]}$ and $\sqrt{V a r [\hat{α}]}$ , biases $(E [\hat{k}] - k)$ and $(E [\hat{α}] - α)$ .

Recently, an ANN estimator based on an ANN has been proposed, which inherits some advantages of RSK and XU estimators by using five feature parameters from ultrasound backscatter envelope signals as its input. As shown in Zhou et al.,³⁶ compared to the RSK estimator and XU estimator, the ANN estimator performs better in terms of NSD and RRMSE of the parameters to be estimated and the calculation speed. However, the disadvantage of an ANN is that its accuracy of parameter estimation is greatly affected by the sample size N. Within a certain range, the larger the size of the samples, the better is the performance of parameter estimation. As the ultrasonic backscatter envelope signal is generally not stationary and ergodic, its statistical characteristics may vary with the actual measurement position. A larger sample size corresponds to a larger actual measurement size, so using a larger sample size will make quantitative ultrasound images have lower resolution and weaker specificity. In general, we hope that the sample size of each estimation is as small as possible while ensuring that the sample size is statistically significant, and the parameter estimation results have acceptable accuracy.

To evaluate the impact of the sample size on the ANN estimator, we followed the method in Zhou et al.,³⁶ to estimate the parameter of the HK distribution when the sample sizes were 100, 1000, 5000, and 10,000 in the domains of $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ . The training sets of i.i.d. envelope signal samples of the HK distribution were generated for $k \in {0.01, 0.02, \dots, 0.99, 1.00}$ and $α \in {0.1, 0.12, \dots, 9.9, 9.0}$ , which were used as target sets. The sets were generated for $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ , which were used as reference standards for the test sets. The ANN architecture is shown in Figure 1, which is a four-layer feed-forward backpropagation (BP) network. Furthermore, we explored the performance of the proposed MLE1 for parameter estimation when the sample size was small (i.e., N = 100) in the domains of $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ . Here, the numerical integration method based on CAPSO was used to calculate the PDF of the HK distribution. In the pattern search algorithm, the initial values of the parameters ɛ $and α$ which are estimated, were set to 4 and 10, respectively. Their lower boundaries were set to 0.01 and 0, and their upper boundaries were set to 4 and 10, respectively.

Estimation Based on Simulated Ultrasound Images

Experimental setup

To verify the performance of the proposed MLE1 in the ultrasonic envelope signal, the ultrasound simulation software Field II was used to simulate a series of ultrasound images. Field II software is an ultrasound simulation software based on the MATLAB platform, which generates ultrasound echo signals from the scatterer phantoms using the linear acoustics theory. A linear array transducer with a focal length of 20 mm was used, and it was excited with a Hanning windowed sinusoidal pulse. The number of physical elements was 512 and 64 for the active elements. The width and height of the element were 0.77 and 5 mm, respectively. Its center and sampling frequency frequencies were 10 and 100 MHz, respectively. The sound speed and wavelength were 1540 m/s and $1.54 \times 10^{- 4}$ m, and the number of lines in the images was 50. The unfiltered and uncompressed RF data was then decimated along the beam axis by a factor of 10 to obtain i.i.d. samples modeled by a single HK distribution. The ultrasound echo envelope was computed as the norm of the Hilbert transform of the RF data.

We generated a three-dimensional computational phantom with dimensions of 20 mm (axial), 20 mm (lateral), and 2 mm (elevation). The center of the volume was located at the geometric focus of the transducer. The resolution cell was considered to have an ellipsoidal shape with diameters equal to the resolutions in each dimension: 0.154 mm (axial), 0.64 mm (lateral), and 0.64 mm (elevation) which was obtained by measuring the width of the envelope of the Point-Spread Function at −6 dB (Full Width at Half Maximum). This configuration leads to a resolution cell with a volume of 0.2642 mm³, based on an ellipsoid model. Thus, within the computational phantom, the average number of scatterers per unit volume (i.e., the density of scatterers) of N_s per resolution cell can be formulated as $N_{s} \times 3028$ scatterers within the numerical phantom.

Generation of simulation phantom and parameters estimation

To perform the simulations, a set of phantoms with various scatterer densities was scanned to generate the simulated ultrasonic echo signals. The signals containing only the diffuse component were simulated by placing randomly located scatterers in the phantom volume at spatial positions distributed according to a uniform distribution. The amplitudes of all scatterers were distributed according to a normal distribution of mean 0 and variance 1. Theoretically, the coherent component is related to the strength of the specular reflection or the periodic organization of the scatterers. Hence, periodically located scatterers along the transducer axis were used to simulate the ultrasonic echo signal containing a coherent component. The coordinates in the plane perpendicular to the transducer axis were randomly placed according to a uniform distribution, whereas their coordinates along that axis were separated by a constant distance, $λ / 2$ , where $λ$ denotes the wavelength of the ultrasound. The amplitude of the periodic scatterers was fixed to a constant of 1 to match the average intensity value, $E [A^{2}]$ , of the randomly located scatterers.^32,35

Sequences of ultrasound images were simulated by randomly located scatterer phantoms with the number of scatterers per resolution cell varying from 1 to 10, which only contain diffuse components. A total of 60 simulated ultrasonic images were generated for each scatterer density value. To simulate phantoms that can generate ultrasonic echo signals that contain not only diffuse components but also coherent components, a fixed density of five randomly located scatterers per resolution cell was considered, and the scatterer density for periodic scatterers was varied from 1 to 5 scatterers per resolution cell with a step of 1. A total of 60 images were simulated for each scatterer density value.

For each simulated sequence, the simulated ultrasonic envelope signals were estimated using the MLE1, RSK, XU, and ANN. The estimated parameter values of α and 1/(k + 1) were averaged over 60 images within a region of interest (ROI) . The ROI was a rectangular area centered at the geometric focus of the simulation phantom, and the size of sample N was 2000. The parameter α is called the effective number of scatterers and can be expressed as $α_{s} N_{*}$ , where $N_{*}$ is the actual number of randomly located scatterers per resolution cell and α_s is a parameter describing the lack of uniformity of the scattering cross sections in the range cell. Moreover, extremely small values of α_s correspond to strong fluctuations in the scattering cross section, and large values of α_s correspond to weak fluctuations. The value of α_s varies from 0 to 1 for the point scatterer phantom. Assuming a K distribution, we can obtain the value of α_s using the methods described in Blacknell and Tough⁴⁴ and Narayanan et al.⁴⁵ The value of constant α_s was .94, which was obtained by averaging all α_s estimated from the ROI of the simulated ultrasonic images. For the pattern search algorithm of the proposed MLE1, in the estimation of the ultrasound envelope signal containing only the diffuse component, the initial values of the parameters σ, ɛ, $and$ α, which are estimated, were set to 5, 0, and 10, respectively. Their lower boundaries were set to 0, 0, and 0, and their upper boundaries were set to 5, 1, and 10, respectively. In the estimation of the ultrasound envelope signal containing both diffuse and coherent components, however, the initial values of the parameters σ, $ε, and α$ which are estimated, were set to 5, 10, and 10, respectively. Their lower boundaries were set to 0, 0, and 0, and their upper boundaries were set to 5, 15, and 10, respectively. Here, the numerical integration method based on CAPSO was used to calculate the PDF of the HK distribution.

Results and Discussion

Computer Simulation of Independent and Identically Distributed Samples From HK Distributions

The total absolute RRMSE, RB, and NSD obtained from the parameters that were estimated by MI, MA, RSK, XU, and the proposed MLE1 in the domains of $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ are shown in Table 2, based on 1000 estimates. Note that the results of MI, MA, RSK, and XU are quoted from those in Destrempes et al.³⁵ Moreover, Figures 6 and 7 depict the plots of RRMSE, RB, and NSD between the parameters estimated by the proposed MLE1 and the reference standard parameters for N = 1000, based on 1000 computer simulation estimates. We can see from Table 2 that the proposed MLE1 is the best in terms of RRMSE, RB, and NSD, especially for the estimation of parameter α. In addition, compared to using the simplest traditional trapezoidal integration method to calculate numerical integration in the proposed MLE1, we find that the use of the advanced numerical integration method based on CAPSO has almost no obvious advantage. Specifically, the MLE1 which adopts the numerical integration method based on CAPSO has better estimation of parameter k, but worse estimation of parameter α. This result further supports the point made in Section 2.5. Considering that the estimation accuracy of the proposed MLE1 for parameter α is significantly better than that for parameter k on the whole, to obtain more accurate parameter estimation for parameter k, the numerical integration method based on CAPSO is adopted to calculate the PDF of HK distribution in subsequent simulation experiments. Of course, one can consider the traditional trapezoidal numerical integration method in terms of computational costs. Furthermore, as shown in Figures 6 and 7, a relatively large estimation error of parameters α was observed in the area of the parameter domain where the value of the reference parameter α was small. In addition, when the value of the reference parameter k was low, it was difficult to obtain an accurate estimate of parameter k.

Table 1.

Main Parameters and Statistics Discussed in This Paper.

Definition	Notion
α	Scatter clustering parameter
ɛ ²	Coherent signal power
$2 σ^{2} α$	Diffuse signal power
$k = \sqrt{2 κ}$	Ratio of the coherent to diffuse signal
$\bar{I}$	Mean intensity
A	Amplitude of envelope signal
$κ$	Structure parameter
$X_{H K}$	X-statistics
$U_{H K}$	U-statistics
$R_{v}$	Signal-to-noise ratio (SNR)
$S_{v}$	Skewness
$K_{v}$ K v	Kurtosis
$L (\bar{θ})$	Log-likelihood function

Table 2.

Improvements in Bias and Variance of Estimators. The Sample Size is N = 1000. The Parameter Domain is $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ . Estimation Methods: (1) MI; (2) MA; (3) XU; (4) RSK; (5) MLE1₁ (Using Numerical Integration Method Based on CAPSO); (6) MLE1₂ (Using Trapezoidal Numerical Integration Method).

	MI	MA	XU	RSK	MLE1₁	MLE1₂
Total absolute value of the relative bias of α	47.9	19.6	17.5	23.9	13.2	4.7
Total absolute value of the relative bias of k	42.2	33.2	34.7	60.7	24.3	28.0
Total normalized standard deviation of α	108.9	73.7	70.3	84.5	8.2	8.3
Total normalized standard deviation of k	117.0	92.0	91.1	62.1	32.6	48.7
Total value of the RRMSE of α	120.3	76.7	72.7	88.2	16.8	10.4
Total value of the RRMSE of k	126.6	99.1	99.1	92.7	45.3	59.6

Figure 8 shows the biases and standard deviations for the estimated parameters using the proposed MLE1 in the domains of $k \in {0.0, 0.1, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 79, 80} .$ The sample size was N = 1000, and each value was obtained by averaging 100 estimations. Overall, we found that the proposed MLE1 possesses the capability to estimate the parameters from sparse and concentrated media. Compared with the results of XU estimator in Cristea et al.⁷ (i.e., Figures 2 and 3), the proposed MLE1 had smaller biases and standard deviations in estimating the parameters of media with more than 10 scatterers per resolution cell (i.e., the average biases and average variance were 2.156 and 3.560 for parameter α, and are 0.153 and 0.188 for parameter k, respectively). A positive bias was observed in the parameter domain where the value of α was small and k was large, whereas in other parameter domains, there was a negative bias that became stronger as α increased and k decreased. The variation trend of the bias of parameter α in the parameter domain was consistent with that obtained by the XU estimator in Cristea et al.⁷ (Figure 3). However, the variation range of the bias of the proposed MLE1 was smaller, especially when the value of α was large and the sample size N was 1000. The range of the bias of α calculated by the MLE1 was between approximately −2 and −8. When calculated by the XU estimator, the bias range of α was between approximately −20 and −50. For the estimation of parameter k, when the value of k was larger and the value of α was smaller, it had a strong negative bias. In general, as α increases, the bias of k gradually varies from a negative value to a positive value. Additionally, the standard deviation of α estimated by the proposed MLE1 from simulated i.i.d. samples increased with α, while there was no obvious trend of change for the standard deviation of k. Moreover, Cristea et al.⁷ shows that the biases and standard deviations of the parameter α estimated by the XU estimator decreased with an increase in the sample size. However, even when the sample size was 50,000, the result of using the XU estimator was still worse than that of using the proposed MLE1 when the sample size was 1000.

Table 3 shows the mean RB, NSD, and RRMSE of the parameter estimated by the ANN estimator in the domains of $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ for the four groups of sample sizes (N = 100, 1000, 5000, 10,000), based on 100 estimates of computer simulations. We found that the performance of parameter estimation using an ANN estimator greatly depended on the sample size; that is, the larger the sample size, the more accurate the parameter estimation. This variation trend of the results agrees with that in Zhou et al.³⁶ (Table 2–4). The RB, NSD, and RRMSE of the proposed MLE1 are presented in Figure 9 for the domains $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ when the sample size N is 100. Using the proposed MLE1, the variation trends of RB, RRMSE, and NSD parameters α and k of the sample size of 100 was similar to that of the sample size of 1000. In addition, comparing the results of Tables 2 and 3, we found that when the sample size was reduced 10 times, the performance of parameter estimation using the proposed MLE1 was only slightly worse. This shows that the proposed MLE1 can effectively estimate the parameters even when the sample size is small, which has potential advantages in quantitative ultrasound.

Table 3.

Total Absolute Value of RB, NSD, and RRMSE Between the Parameter Estimated by the ANN Estimator or MLE1, and the Reference Standard of α and k for Different Sample Size, N, Based on 100 Estimates of Computer Simulations. The Parameter Domain is $k \in {0.1, 0.2, \dots, 0.9, 1.0}$ and $α \in {1, 2, \dots, 9, 10}$ .

Estimators	ANN	ANN	ANN	ANN	MLE1
Sample size N	100	1000	5000	10000	100
Total absolute value of the relative bias of α	58.1	25.4	16.0	12.5	14.5
Total absolute value of the relative bias of k	79.6	74.7	73.6	70.9	28.6
Total normalized standard deviation of α	36.5	22.4	15.6	10.5	16.7
Total normalized standard deviation of k	13.8	13.2	14.2	12.9	48.4
Total value of the RRMSE of α	66.5	30.9	19.0	14.9	23.9
Total value of the RRMSE of k	81.7	76.5	73.5	72.9	61.3

Estimation Based on Simulated Ultrasound Images

Examples of simulated B-mode images corresponding to different concentrated simulated phantoms, described in Section 3.3, are illustrated in Figure 10, where log-compression and downsampling were applied to the simulated ultrasonic echo signal for visualization. Figure 11 shows the variation curve of the theoretical parameter values and the corresponding estimated parameter values obtained from the ultrasonic echo signal containing only the diffuse component using different estimators. As can be seen from this figure, the scatterer clustering parameter α increased with the growth of the average number of randomly located scatterers per resolution cell, and the results of parameter estimation using MLE1 are closer to corresponding theoretical values than that using other estimators. Furthermore, examples of the histogram of the ultrasonic echo envelope signals that contain both diffuse and coherent components and the probability density function curve of the HK distribution plotted by the parameters estimated by using the MLE1 from the corresponding envelope signals are shown in Figure 12. It can be observed that the curve drawn by the estimated parameters fits the histogram envelope well. The mean RB, NSD, and RRMSE for α were .1176, .0596, and .1357, respectively, and for k were .2847, .1686, and .3678, respectively.

Note that the theoretical value of k is calculated as $k = \sqrt{2 κ}$ , where $κ$ is equal to $\frac{ε^{2}}{2 α σ^{2}}$ and can be regarded as the ratio of periodic to randomly located scatterers in this simulation experiment. In conclusion, the aforementioned simulation experiment results show that the proposed MLE1 can reliably estimate simulated ultrasound images.

Discussion

In this paper, we systematically investigated the application of MLE in parameter estimation of HK distribution and proposed a new MLE for the HK distribution on the basis of overcoming the mistakes or deficiencies in the conventional MLE. We conclude that it is more reasonable to use equation (4) to calculate the PDF of HK distribution, and the performance of numerical integration method has no obvious influence on the accuracy of parameter estimation. One can choose specific numerical integration methods based on actual demand. Moreover, the results obtained from the simulation experiments using simulated i.i.d. samples from the HK distribution and simulated ultrasonic echo signals confirmed the feasibility of using the proposed MLE1 to obtain accurate parameter estimation for the HK distribution. In addition, compared to using the XU estimator, the proposed MLE1 had a smaller bias and variance in the parameter estimation of concentrated media. This advantage will greatly contribute to the diagnosis of diseases when the tissue has a high number of scatterers per resolution cell based on statistical QUS technology, such as the classification of cancerous lymph nodes. Furthermore, there was a weak dependency between the performance of the proposed MLE1 and the sample size, and a relatively accurate parameter estimation can be obtained even when the sample size is small. This is a considerable advantage of the proposed MLE1 over the latest XU estimator and ANN estimator because the hypothesis that speckles over a wide area of tissues scanned in clinical disease detection that have the same statistical properties may fail. Although there are several methods to increase the number of uncorrelated samples, as mentioned in Cristea et al.,⁷ the parameter estimation method based on a small sample size is more popular because it can effectively improve the resolution and specificity of quantitative ultrasound images and thus enhance the accuracy of disease diagnosis. Moreover, as shown in Figure 11, in the simulation experiments based on ultrasonic echo signals, the proposed MLE1 yields a more reliable and stable HK parameter estimation compared to other estimators. In conclusion, compared with other estimators, the MLE1 proposed in this paper is more accurate and versatile. It can characterize in a wide range of number density of scatterers and is less dependent on sample size. Therefore, it will further improve the accuracy and applicability of the statistical QUS technology based on HK distribution.

It is well-known that MLE has the disadvantage of being more time-consuming than the methods of moments and ANNs. In the case of the first experiment detailed in Section 3.2, one estimation took an average of 4s (using the traditional trapezoidal numerical integration method) and 8s (using the numerical integration method based on CAPSO). The computational load of the MLE1 is determined by the selection of the numerical integration method, global optimization algorithm, and sample size. In addition, the accuracy of parameter estimation using the proposed MLE1 is significantly affected by the global optimization algorithm. It is worth mentioning that the selection of the initial values and upper and lower boundaries of the parameters to be estimated in the pattern search global optimization algorithm used in this study plays an important role in the performance of the MLE1. Therefore, we can use some HK distribution parameter estimators with fast operation speed, such as XU and ANN estimators, to first estimate the HK distribution parameters, and then set the estimated parameter value as the initial value, and set the upper and lower boundaries of the parameters to be estimated near the initial value. This operation may make the MLE1 faster and more accurate. Moreover, the pattern search algorithm used in this study was determined by comparing several common gradient-free global optimization algorithms, such as the particle swarm optimization algorithm and genetic algorithm, under a trade-off between the calculation cost and parameter estimation accuracy. In future work, new and effective numerical integration methods and global optimization algorithms can be used to increase the accuracy and efficiency of the MLE1.

Note that the scatterers used in the field II ultrasound simulation are treated as dimensionless points with a position and amplitude. This hypothesis is inappropriate for biological tissues, as described in Section 5.4 of Destrempes et al.³⁵ Furthermore, the validation of the proposed MLE1 in this study is based on simulation data, which is a limitation of our study. Therefore, in future work, more clinical experiments can be conducted to further validate the performance of the proposed MLE1. Moreover, this paper only discusses the usefulness of the proposed MLE1 and XU estimator for the characterization of both sparse and concentrated media, more comprehensive comparative experiments involving more estimators should be performed to further evaluate the performance of each estimator.

Conclusions

In this study, we proposed an effective parameter estimation method of the HK distribution based on the MLE after analyzing the deficiencies of conventional MLE of HK distribution. The simulation experiment results confirmed the feasibility of the proposed MLE1 and indicated that the proposed MLE1 is more accurate and robust than other estimators when estimating parameters α and k of the HK distribution. Thus, we can conclude that the proposed MLE1 has potential practical value for ultrasonic parameter imaging in tissue characterization.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Grants from the Yunnan Fundamental Research Projects (No. 202101AS070031), and University Key Lab of Electronic Information Processing of High Altitude Medicine, Yunnan Province.

ORCID iD

Yufeng Zhang

References

Byra

Nowicki

Wróblewska-Piotrzkowska

Dobruch-Sobczak

Classification of breast lesions using segmented quantitative ultrasound maps of homodyned K distribution parameters. Med Phys. 2016;43(10):5561-9.

Zhou

Zhang

Tsui

PH.

Hepatic steatosis assessment using quantitative ultrasound parametric imaging based on backscatter envelope statistics. Appl Sci. 2019;9(4):661.

Zhou

Gao

Zhang

Tsui

PH.

Ultrasound backscatter envelope statistics parametric imaging for liver fibrosis characterization: a review. Ultrason Imaging. 2020;42(2):92-109.

Mamou

Oelze

ML.

Quantitative Ultrasound in Soft Tissues. Dordrecht: Springer Netherlands; 2013, pp. 1-25, 219-52.

Zhou

Jia

Tsui

PH.

A review of ultrasound tissue characterization with mean scatterer spacing. Ultrason Imaging. 2017;39(5):263-82.

Prager

Gee

Treece

Berman

LH.

Decompression and speckle detection for ultrasound images using the homodyned K-distribution. Pattern Recognit Lett. 2003;24(4-5):705-13.

Cristea

Collier

Franceschini

Mamou

Cachard

Basset

Quantitative assessment of media concentration using the homodyned K distribution. Ultrasonics. 2020;101:105986.

Lizzi

Ostromogilsky

Feleppa

Rorke

Yaremko

MM.

Relationship of ultrasonic spectral parameters to features of tissue microstructure. IEEE Trans Ultrason Ferroelectr Freq Control. 1987;34(3):319-29.

Lizzi

Alam

Mikaelian

Lee

Feleppa

EJ.

On the statistics of ultrasonic spectral parameters. Ultrasound Med Biol. 2006;32(11):1671-85.

10.

Oelze

O’Brien

Jr Blue

Zachary

JF.

Differentiation and characterization of rat mammary fibroadenomas and 4T1 mouse carcinomas using quantitative ultrasound imaging. IEEE Trans Med Imaging. 2004;23(6):764-71.

11.

Oelze

Mamou

Review of quantitative ultrasound: envelope statistics and backscatter coefficient imaging and contributions to diagnostic ultrasound. IEEE Trans Ultrason Ferroelectr Freq Control. 2016;63(2):336-51.

12.

Lizzi

Greenebaum

Feleppa

Elbaum

Coleman

DJ.

Theoretical framework for spectrum analysis in ultrasonic tissue characterization. J Acoust Soc Am. 1983;73(4):1366-73.

13.

Tsui

Yeh

Liao

Chang

Kuo

Chang

, et al. Ultrasonic Nakagami imaging: a strategy to visualize the scatterer properties of benign and malignant breast tumors. Ultrasound Med Biol. 2010;36(2):209-17.

14.

Hao

Bruce

Pislaru

Greenleaf

JF.

Characterization of reperfused infarcted myocardium from high-frequency intracardiac ultrasound imaging using homodyned K distribution. IEEE Trans Ultrason Ferroelectr Freq Control. 2002;49(11):1530-42.

15.

Gennisson

Deffieux

Fink

Tanter

Ultrasound elastography: principles and techniques. Diagn Interven Imaging. 2013;94(5):487-95.

16.

Urban

Chen

Fatemi

A review of shearwave dispersion ultrasound vibrometry (SDUV) and its applications. Curr Med Imaging Rev. 2012;8(1):27-36.

17.

Anavekar

JK.

Doppler echocardiography: a contemporary review. J Cardiol. 2009;54(3):347-58.

18.

Harris

Kuppurao

Quantitative Doppler echocardiography. BJA Educ. 2016;16(2):46-52.

19.

Lee

Han

JK.

Evaluation of hepatic steatosis by using acoustic structure quantification US in a rat model: comparison with pathologic examination and MR spectroscopy. Radiology. 2017;285(2):445-53.

20.

Zhou

Lin

Fang

Liu

Tsui

PH.

Three-dimensional visualization of ultrasound backscatter statistics by window-modulated compounding Nakagami imaging. Ultrason Imaging. 2018;40(3):171-89.

21.

Fang

Zhou

Chang

Wan

Tsui

PH.

Ultrasound parametric imaging of hepatic steatosis using the homodyned-K distribution: an animal study. Ultrasonics. 2018;87:91-102.

22.

Zhou

Wan

Tsui

PH.

A computer-aided diagnosis scheme for detection of fatty liver in vivo based on ultrasound kurtosis imaging. J Med Syst. 2016;40(1):33-9.

23.

Zhou

Tai

Wan

Tseng

Lin

, et al. Hepatic steatosis assessment with ultrasound small-window entropy imaging. Ultrasound Med Biol. 2018;44(7):1327-40.

24.

Destrempes

Cloutier

A critical review and uniformized representation of statistical distributions modeling the ultrasound echo envelope. Ultrasound Med Biol. 2010;36(7):1037-51.

25.

Shankar

PM.

A model for ultrasonic scattering from tissues based on the K distribution. Phys Med Biol. 1995;40(10):1633-49.

26.

Mohana Shankar

. A general statistical model for ultrasonic backscattering from tissues. IEEE Trans Ultrason Ferroelectr Freq Control. 2000;47(3):727-36.

27.

Jakeman

Tough

RJA

. Generalized K distribution: a statistical model for weak scattering. J Opt Soc Am A. 1987;4(9):1764-72.

28.

Tsui

Wan

Tai

Shu

YC.

Effects of estimators on ultrasound Nakagami imaging in visualizing the change in the backscattered statistics from a Rayleigh distribution to a pre-Rayleigh distribution. Ultrasound Med Biol. 2015;41(8):2240-51.

29.

Trop

Destrempes

El Khoury

Robidoux

Gaboury

Allard

, et al. The added value of statistical modeling of backscatter properties in the management of breast lesions at US. Radiology. 2015;275(3):666-74.

30.

Destrempes

Trop

Allard

Chayer

Garcia-Duitama

El Khoury

, et al. Added value of quantitative ultrasound and machine learning in BI-RADS 4–5 assessment of solid breast lesions. Ultrasound Med Biol. 2020;46(2):436-44.

31.

Byra

Kruglenko

Gambin

Nowicki

Temperature monitoring during focused ultrasound treatment by means of the homodyned K distribution. Acta Phys Pol A. 2017;131(6):1525-8.

32.

Dutt

Greenleaf

JF.

Ultrasound echo envelope analysis using a homodyned K distribution signal model. Ultrason Imaging. 1994;16(4):265-87.

33.

Martin-Fernandez

Cardenes

Alberola-Lopez

. Parameter estimation of the homodyned K distribution based on signal to noise ratio. In: IEEE Ultrasonics Symposium, pp. 981-4. doi:10.1016/j.ultras.2020.106308.

34.

Hruska

Oelze

ML.

Improved parameter estimates based on the homodyned K distribution. IEEE Trans Ultrason Ferroelectr Freq Control. 2009;56(11):2471-81.

35.

Destrempes

Porée

Cloutier

Estimation method of the homodyned K-distribution based on the mean intensity and two log-moments. SIAM J Imaging Sci. 2013;6(3):1499-530.

36.

Zhou

Gao

Tai

Tseng

, et al. Parameter estimation of the homodyned K distribution based on an artificial neural network for ultrasound tissue characterization. Ultrasonics. 2021;111:106308.

37.

Lin

Cheng

Huang

Lin

Wan

Tsui

PH.

Detecting changes in ultrasound backscattered statistics by using Nakagami parameters: comparisons of moment-based and maximum likelihood estimators. Ultrasonics. 2017;77:133-43.

38.

Redner

Walker

HF.

Mixture densities, maximum likelihood and the EM algorithm. SIAM Rev. 1984;26(2):195-239.

39.

Zhang

Deng

Cai

Zhang

Zhou

, et al. Assessment of homodyned K distribution modeling ultrasonic speckles from scatterers with varying spatial organizations. J Healthc Eng. 2017;2017:8154780.

40.

Wang

Zhang

Zhao

, et al. A numerical integral algorithm based on the CAPSO to improve the estimation for the parameters of the homodyned-K distribution. In: 2018 IEEE 18th International Conference on Communication Technology (ICCT), pp. 1342-45. doi:10.1109/ICCT.2018.8599949.

41.

Natick. Global Optimization Toolbox User’s Guide (R2020b). The MathWork Inc. 2-3:155-250.

42.

Zhu

. An adaptive particle swarm optimization algorithm based on Cloud Model. Adv Mater Res. 2010;129-131:612-6.

43.

Solving numerical integration by particle swarm optimization. In: Zhu

Zhang

Liu

, eds. International Conference on Information Computing and Applications. Berlin, Heidelberg: Springer; 2010, pp. 228-35.

44.

Blacknell

Tough

RJ.

Parameter estimation for the K-distribution based on [z log(z)]. IEE Proc Radar Sonar Navig. 2001;148(6):309-12.

45.

Narayanan

Shankar

Reid

JM.

Non-Rayleigh statistics of ultrasonic backscattered signals. IEEE Trans Ultrason Ferroelectr Freq Control. 1994;41(6):845-52.