Anca Cristea

Nicolas Collier

Emilie Franceschini

Jonathan Mamou

Christian Cachard

Olivier Basset

Quantitative assessment of media concentration using the Homodyned K distribution

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

resolution cell and/or highly variable scattering cross-sections is the K distribution [START_REF] Jakeman | On the statistics of K-distributed Noise[END_REF], whereas the Homodyned K distribution [START_REF] Hao | Characterization of reperfused infarcted myocardium from high-frequency intracardiac ultrasound imaging using Homodyned K distribution[END_REF], [START_REF] Trop | The added value of statistical modeling of backscatter properties in the management of breast lesions at US[END_REF], [START_REF] Jakeman | Generalized K distribution: a statistical model for weak scattering[END_REF] covers all the configurations, including the presence of a coherent component, and is thus the most general model used in ultrasound. A simpler model is the Nakagami distribution [START_REF] Tsui | Ultrasonic Nakagami imaging: A strategy to visualize the scatterer properties of benign and malignant breast tumors[END_REF], [START_REF] Shankar | A general statistical model for ultrasonic backscattering from tissues[END_REF], [START_REF] Shankar | Classification of ultrasonic B-mode images of breast masses using Nakagami distribution[END_REF], [START_REF] Hampshire | Probability density of myocardial ultrasonic backscatter[END_REF], [START_REF] Nakagami | Statistical characteristics of short-wave fading[END_REF], [START_REF] Nakagami | The m distribution-a general formula of intensity distribution in rapid fading[END_REF], which has been used extensively in tissue characterization applications because of the robustness and simplicity of the associated parameter estimation methods. Comparative studies using the aforementioned models, as well as less common ones: Weibull, Generalized Gamma, Lognormal or other generalizations have been conducted on ultrasound data acquired from the human heart [START_REF] Eltoft | Modeling the Amplitude Statistics of Ultrasonic Images[END_REF], human skin [START_REF] Raju | Statistics of envelope of high-frequency ultrasonic backscatter from human skin in vivo[END_REF], [START_REF] Pereyra | Modeling ultrasound echoes in skin tissues using symmetric α-stable processes[END_REF] and human lymph nodes [START_REF] Mamou | Three dimensional high-frequency backscatter and envelope quantification of cancerous human lymph nodes[END_REF], [START_REF] Bui | Modeling the envelope statistics of three-dimensional high-frequency ultrasound echo signals from dissected human lymph nodes[END_REF]. A rigorous and comprehensive review of the theoretical models has been published in [START_REF] Destrempes | A critical review and uniformized representation of statistical distributions modeling the ultrasound echo envelope[END_REF].

The Homodyned K and the Nakagami distributions are the most frequently used models for ultrasound signal envelopes. Both distributions are sensitive to the number of scatterers per resolution cell through their shape parameters (the scatterer clustering parameter α for the Homodyned K and the m parameter for the Nakagami). Moreover, the Nakagami distribution can be viewed as a good approximation of the Homodyned K distribution for a range of parameter values compatible with those measured in some tissue characterization applications [START_REF] Destrempes | Review of Envelope Statistics Models for Quantitative Ultrasound Imaging and Tissue Characterization[END_REF]. In recent years, comparative studies have presented the Homodyned-K parameter α as more appropriate than the Nakagami parameter m for describing media with high numbers of scatterers per resolution cell, such as densely packed cells in tumors [START_REF] Oelze | Statistics of scatterer property estimates[END_REF], and as the best statistical QUS estimate for the classification of cancerous lymph nodes [START_REF] Mamou | Three dimensional high-frequency backscatter and envelope quantification of cancerous human lymph nodes[END_REF]. However, its use has been more limited because of the analytical complexity of the model and the difficulty of obtaining reliable estimates in situations where the number of scatterers per resolution cell is superior to 10 [START_REF] Oelze | Statistics of scatterer property estimates[END_REF].

The purpose of the present paper is twofold. Firstly, we assess the reliability of the Homodyned K scatterer clustering parameter α estimates computed using the recently published XU estimator [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF], for a range of values compatible with measurements obtained from both sparse and concentrated media. We consider "concentrated media" to be media which contains more than 10 scatterers per resolution cell when imaged at a specific resolution. Secondly, we assess the ability of the parameter α to evaluate scattering properties from sparse and concentrated media, based on computer simulations and in vitro experiments. Previous simulation studies [START_REF] Dutt | Ultrasound echo envelope analysis using a Homodyned K distribution signal model[END_REF], [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF], [START_REF] Hruska | Improved parameter estimates based on the Homodyned K distribution[END_REF] have assumed dimensionless scatterers and in general focused on low numbers of scatterers per resolution cell ( [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF], [START_REF] Hruska | Improved parameter estimates based on the Homodyned K distribution[END_REF]), thus limiting the interpretation of the results to the case of a low density of uniformly distributed scatterers of sufficiently small size (i.e., sparse media). The authors of [START_REF] Dutt | Ultrasound echo envelope analysis using a Homodyned K distribution signal model[END_REF] study concentrations up to 35 scatterers/resolution cell, but do not discuss the variance of the estimates, thus offering an additional motivation for our study. In our study, we use ultrasound simulation software [START_REF] Varray | CREANUIS: A non-linear radiofrequency ultrasound image simulator[END_REF] to calculate the interaction between a predefined ultrasound pulse and distributions of non-overlapping scatterers of a given size, for a more realistic recreation of experimental conditions, including the consideration of correlations that arise between scatterer positions at high concentrations (up to 68 scatterers per resolution cell). Regarding the experimental use of the XU estimator, the only in-vitro experiment performed on concentrated media was conducted on aggregating red blood cells [START_REF] Destrempes | Unifying concepts of statistical and spectral quantitative ultrasound techniques[END_REF], with the purpose of demonstrating the relationship between the coherent and diffuse signal power determined via Homodyned-K modeled statistics, on one hand, and BSC-related parameters, on the other hand. Our study focuses on the relationship between the parameter α and the number of scatterers in tissue-mimicking phantoms, of which the concentrated phantoms (up to 101 scatterers per resolution cell) are used to mimic densely packed cells in tumors [START_REF] Franceschini | Experimental assessment of four ultrasound scattering models for characterizing concentrated tissue-mimicking phantoms[END_REF]. The diameters of the phantom scatterers and the simulated scatterers are identical, in order to facilitate comparisons. In addition, the analogous Nakagami parameter m is also estimated throughout the study, for comparative purposes.

We present the models and methods for both simulation and experiments in detail in Section II. Section III highlights the results, followed by a discussion in Section IV and conclusions in Section V.

II. MODELS AND METHODS

A. Models: The Homodyned K and Nakagami Distributions

The Homodyned K model is defined by the following probability density function [START_REF] Destrempes | A critical review and uniformized representation of statistical distributions modeling the ultrasound echo envelope[END_REF]:

( ) ( ) ( ) 2 2 0 0 0 1 . 2 HK x A A xJ x Ax P J dx α σ ε - ∞   = +     ∫ (1)
denotes the Bessel function of first kind and order 0, α is the scatterer clustering parameter, 2σ 2 α is the diffuse signal power and ε 2 is the coherent signal power. The analytical complexity of the distribution and the absence of a closed-form expression make parameter estimation particularly difficult. Moment-based methods have been proposed in [START_REF] Dutt | Ultrasound echo envelope analysis using a Homodyned K distribution signal model[END_REF] and [START_REF] Martin-Fernandez | Parameter estimation of the Homodyned K distribution based on signal to noise ratio[END_REF] but their limited precision determined a need to develop more complex methods that rely on other statistics as well as elegant optimization methods. Hruska and Oelze [START_REF] Hruska | Improved parameter estimates based on the Homodyned K distribution[END_REF] proposed a method which uses the values of fractional-order SNRs, skewness and kurtosis to estimate α and the ratio between the coherent and diffuse signal components. The most recently developed method is the "XU estimator" [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF]. Of all methods, the XU estimator delivers the best results in terms of precision and speed, so we choose to use it in our work. We give a brief description of the algorithm here.

The XU estimator uses the mean value of the intensity I (the square of the envelope amplitude) together with two log-moments of the envelope (also called the X and U statistics) as inputs for the estimation. These are defined as:

,

I µ = (2) log log , I I X I I = - (3)
log log ,

U I I = - (4)
where the symbol (•) denotes averaging. In order to develop analytical expressions for the two log-moments as a function of the model parameters, it is convenient to introduce an algorithmic parameter defined as γ=ε 2 /2σ 2 . Consequently, ( , ) HK X α γ and ( , ) HK U α γ ([38] eq. 3.3, eq.3.4]) are obtained. For computational efficiency, the value of α is limited to a maximum value αmax (here set to 80, rounding up the value of 79.5 used in the original work). The estimation is formulated as a minimization problem:

( ) ( ) 2 2
, ,

H K a r g m in U U ε σ α = - (5)
subject to:

, [I] . , HK max E I X X α α = ≤ = (6)
The parameters α and γ which optimally satisfy the equations are estimated using a binary search algorithm, and the parameters ε and σ can then be obtained from the equations: 2 / ( ),

ε µγ γ α = + (7) 2 / (2( ))
.

σ µ γ α = + (8)
Let us now consider the Nakagami model, which is described by the following probability density function [START_REF] Nakagami | The m distribution-a general formula of intensity distribution in rapid fading[END_REF]:

( ) 2 1 2 2 (A) exp . m m Nak m m m p A A m -   = -   Ω Γ Ω   (9)
The two characteristic parameters are m and Ω, and Γ(•) is the gamma function. The shape parameter m (also called Nakagami parameter) increases with the number of scatterers per resolution cell up to a value of around 1 for fully developed speckle, provided that the coherent component vanishes. The scale parameter Ω is equal to the mean backscattered intensity. Parameter estimation is quite straightforward and is often performed using expressions derived from the method of moments:

2 2 2 [ ] , [ ] E A m Var A = (10) 2 [ ]. E A Ω = (11)
In the present study, a maximum likelihood-based estimation method was adopted, because of its increased sensitivity to small values of the parameter m. The likelihood equation for m:

2 1 1 ln( ) ln( ) ln( ) ( ) 0, N j j m X m N ψ = -Ω + - = ∑ (12)
where ψ(•) represents the digamma function and N the number of samples, is solved by numerical optimization, using the Newton-Raphson method. The initial solution is obtained using [START_REF] Brand | Monitoring of cell death in epithelial cells using high frequency ultrasound spectroscopy[END_REF], and the final solution is reached when the difference between consecutive estimates becomes inferior to 10 -8 . Note that the maximum likelihood estimate of Ω is equivalent to [START_REF] Meziri | In vitro chronic hepatic disease characterization with a multiparametric ultrasonic approach[END_REF].

B. Simulations of Independent and Identically-Distributed (i.i.d.) Samples from Homodyned K Distributions

The first part of the study represents a characterization of the estimation bias and variance of the scatterer clustering parameter α. An extensive exploration of the parameters ε and σ, including our preliminary study presented in [START_REF] Cristea | Quantitative characterization of concentrated cell pellet biophantoms using statistical models for the ultrasound echo envelope[END_REF] has allowed us to conclude that they do not relate exclusively to the number of scatterers per resolution cell and are thus beyond the scope of this paper. We do however mention that the coherent-to-diffuse signal ratio k:

= √ (13)
relates to spatial arrangement of the scatterers and influences the bias and variance of α [START_REF] Dutt | Ultrasound echo envelope analysis using a Homodyned K distribution signal model[END_REF], [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF], [START_REF] Hruska | Improved parameter estimates based on the Homodyned K distribution[END_REF], therefore its role is considered here. High numbers of scatterers per resolution cell (>10) characteristic to concentrated media produce very similar speckle patterns, therefore estimators must be sensitive to small variations in the shape of the distribution in order to extract the correct parameter values. This is a somewhat intuitive explanation for the difficulty of performing precise estimations on concentrated media. A more rigorous explanation for the variance is given by the Cramer-Rao Lower Bound (CRLB) ( [START_REF] Kay | Fundamentals of statistical signal processing: Estimation theory[END_REF], [START_REF] Rice | Mathematical statistics and data analysis[END_REF]):

1 ( ) . ( ) F Var NI α α ≥ ) (14)
The CRLB states that the minimum variance of a non-biased estimate of a distribution parameter (here, α) is inversely proportional to the Fisher information available for the respective parameter and to the number of samples N. The Fisher information describes the sensitivity of the log-likelihood function ln p(A; ) α to the parameter α:

2 2 ln ( ) . F p(A; ) I E α α α   ∂ = -  ∂   (15)
Even if the CRLB does not apply to biased estimates, it does offer a conceptual tool for understanding the evolution of the variance. Note also that the bound does not depend on the estimator, therefore it describes an ideal situation for the given model. As the likelihood function for the Homodyned K distribution is difficult to express in analytical form, we used numerical simulations to observe the evolution of the Fisher information for different values of α and k.

Two sets of simulations of i.i.d. samples were considered in this subsection: one for the Fisher information, and one for the assessment of the estimation bias and variance. The same range of parameters was considered for both sets: α between 1 and 80, and k between 0 and 1 (based on values estimated from ultrasound data). In order to preserve scale, thus removing the influence of the parameter σ, the mean of the intensity was constrained to 1:

+ 2 = 1. ( 16
)
From ( 13) and ( 16), it follows that = /( + 2) and = (1 -)/2 . The i.i.d. samples from a Homodyned K distribution are generated using ( [START_REF] Hruska | Improved parameter estimates based on the Homodyned K distribution[END_REF]):

( ) ( ) 2 2
,

i i i i i x X z Y z ε σ σ = + + (17)
where xi denotes an i.i.d. Homodyned K sample, Xi and Yi are i.i.d. samples of the standard Gaussian distribution and zi is a sample of the Gamma distribution with shape parameter α and scale parameter 1. For the computation of IF(α, k), we averaged 100 results for every (α, k) value pair. Each result was obtained by simulating 10 000 i.i.d. samples and computing the numerical second derivative of the log-likelihood function with respect to α (eq. 15) considering these samples. The likelihoods were evaluated by approximating eq. ( 1) with the power series from [44, eq. 34]. For the assessment of the estimation bias and variance, we performed 100 simulations of 1 000, 10 000 and 50 000 i.i.d. samples for every (α, k) value pair, then performed parameter estimation using the XU estimator. The results are presented in the next section.

C. Simulations of Ultrasound Signals

In order to demonstrate the dependency of the α estimates on the number of scatterers per resolution cell and the number of samples, we have simulated a series of ultrasound signals and performed estimations on collections of these signals.

The first step of the simulation consisted in generating the scatterer positions in a three-dimensional (3D) space. We considered 11 scatterer concentrations between 0.005 and 0.2 in volume fractions (see Table 1). The positions of the scatterer centers were generated randomly from a uniform distribution. As the simulator only considers point scatterers, we enforced that the scatterers cannot overlap by imposing an exclusion distance equal to the scatterer diameter (here set to 13 µm). For each concentration, we generated 10 3D spaces with dimensions: 1 mm (axial) x 5 mm (lateral) x 0.5 mm (elevation), chosen by taking into consideration the imaging frequency, the corresponding wavelength and the duration of the computation.

The second step was the ultrasound signal simulation using the software CREANUIS [START_REF] Varray | CREANUIS: A non-linear radiofrequency ultrasound image simulator[END_REF], designed for efficient nonlinear signal simulation. Even in the absence of signal harmonics, CREANUIS proves to be extremely time-efficient for the simulation of large volumes of data. The simulation method is based on the Born approximation, therefore only considers simple scattering. We considered an array with 64 active elements, a pitch of 68 µm, kerf of 4 µm and height of 2.5 mm, focused at 6 mm. The excitation pulse had a central frequency of 30 MHz and a symmetrical frequency spectrum with a bandwidth measuring 60% of the central frequency. The resolution cell was considered to have an ellipsoidal shape with the diameters equal to the resolutions in each dimension: 0.042 mm (axial), 0.143 mm (lateral) and 0.123 mm (elevation). This configuration results in a resolution cell with a volume of 0.00038 mm 3 . Resolutions were obtained by measuring the width of the envelope of the Point-Spread Function at -6 dB (Full Width at Half Maximum). For each of the 10 volumes, images were taken at 10 locations (equally spaced at 0.05 mm) in the elevation dimension. As the RF signals are spaced at an interval equivalent to the probe pitch, each 1 mm x 5 mm image was composed of 74 RF signals. The sampling frequency was 240 MHz. All images corresponding to a given scatterer concentration were then concatenated along the lateral dimension, resulting in a total number of 7400 RF signals.

In order to have results comparable to those obtained on simulated data, the samples on which the estimation is performed need to be independent. Adjacent lines (in both the lateral and axial dimensions of an US image) are jointly normal distributed, therefore also independent if uncorrelated. Any correlation between adjacent lines can be reduced by downsampling. Here, it was necessary to downsample by a factor of 2 in order to reduce the mean correlation coefficient between adjacent lines to a value of ρ<0.1. The lines were then randomly permuted in order to remove any residual correlation before grouping them for estimation. A downsampling factor of 2 was also applied along the axial dimension, based on the sampling frequency and the maximum frequency of the measured signal spectrum (50 MHz at -40 dB). The average axial correlation was thus also reduced to ρ<0.1. The total equivalent number of independent samples thus obtained for each of the concentrations is approximately 500 000. Finally, the 500 000 samples were divided into groups of different sizes, and we obtain 500 estimates on 1 000 i.i.d. samples, 50 estimates on 10 000 i.i.d. samples and 10 estimates on 50 000 i.i.d. samples.

Simulations

D. Acquisitions from tissue-mimicking Phantoms

The tissue-mimicking phantoms consisted of polyamide microspheres with a mean diameter of 13 µm (Orgasol 2001 EXD NAT1, Arkema, France) gently stirred in a mixture of water and surfactant (Tergitol NP-40, Sigma-Aldrich, France) to prevent sedimentation. The addition of surfactant (less than 0.5%) allows us to avoid the formation of particle clusters. All phantoms had identical scatterer sizes but different scatterer volume fractions ranging from 0.01 to 0.25. The polyamide microspheres have known acoustical properties (sound speed of 2300 m/s, density of 1030 kg/m 3 and Poisson's ratio of 0.42). The concentrated phantoms were designed to mimic the structural properties of densely packed cells in tumors.

As performed previously in [START_REF] Franceschini | Experimental assessment of four ultrasound scattering models for characterizing concentrated tissue-mimicking phantoms[END_REF], a broadband transducer with center frequency 22 MHz, a -6 dB bandwidth of 11-35 MHz, focus of 13 mm and f-number of 2.1 was used for the experiment. The axial, lateral and elevation resolutions are 0.036 mm, 0.140 mm and 0.140 mm, respectively, resulting in a resolution cell volume of 0.00036 mm 3 . The equivalent numbers of scatterers per resolution cell are given in Table 1. The pulse-echo acquisition system was composed of an Olympus model 5073 PR pulserreceiver and a Gagescope model 8500CS acquisition board, and the experimental RF data were registered with a sampling frequency of 250 MHz and 12 bit resolution. The transducer was placed in an agar gel, i.e. a solidified mixture of distilled water and 2% (w/w) agar powder (A9799, Sigma Aldrich, France) so that the distance between the transducer focus and the interface of agar/suspension was equal to 12 mm (see [START_REF] Franceschini | Experimental assessment of four ultrasound scattering models for characterizing concentrated tissue-mimicking phantoms[END_REF]Fig. 4] for a schematic representation of the experiment set-up). The transducer focus was therefore positioned below the agar/suspension interface at a distance of 1 mm. The agar gel helps maintain the focal properties between the transducer and the suspension, by minimizing potential variations in attenuation and speed of sound due to the changes in volume fractions.

The total attenuation coefficients of the tissue-mimicking phantoms were determined using a standard substitution method. The 22 MHz transducer was used in the reflection mode with a reflector on the opposite side. For each phantom, acquisitions of 100 RF signals were performed both with the suspension and with water. The water acquisition was used for normalization. At each acquisition, the total attenuation of the phantoms was calculated using the log spectral difference technique [START_REF] Kuc | Estimating the acoustic attenuation coefficient slope for liver from reflected ultrasound signals[END_REF]. These values were then averaged to obtain the total attenuation coefficients αatt ([dB/MHz/cm]) (Table 1). The high attenuation of the polyamide microspheres causes the signal to decay rapidly within the ROI, thus broadening the envelope distribution and rendering it similar to the distribution of signal envelopes from a sparse medium. This effect, if uncorrected, leads to a drastic underestimation of α, and moreover to a loss of distinction between concentrations. Attenuation compensation is therefore essential.

A total of 6000 RF signals were acquired and stored for the experiment. The region of interest was selected in the focal zone with a rectangular window d1=0.3 mm in length. Point attenuation compensation [START_REF] Oelze | Frequency-dependent attenuation-compensation functions for ultrasonic signals backscattered from random media[END_REF] was applied to the power spectrum of each selected signal in a bandwidth of 11 to 35 MHz, at each depth d and frequency f, assuming the following relationship between the noncompensated and compensated power spectra:

( , ) = ( , )! "# $%% &.(& )* (18)
where the coefficient 8.68 expresses unit conversion from decibel to neper. After attenuation compensation in the frequency domain, the signals are transformed back into the time domain. For simplicity, we assume that the change in frequency content does not determine a significant change in the phase of the signal and thus consider that the phase is not affected by this processing step. The RF signals are not correlated, but downsampling by a factor of 3 is still necessary to reduce the average correlation (to ρ<0.1) in the axial dimension. Approximately 200 000 independent samples are obtained. Three series of estimates are then obtained from groups of samples of different sizes, as it was done in the previous subsection: 200 estimates on 1000 i.i.d. samples, 20 estimates on 10 000 i.i.d. samples and 4 estimates on 50 000 i.i.d. samples.

III. RESULTS

A. Independent and Identically-Distributed Samples from Homodyned K Distributions

Numerical simulations showed that the Fisher information IF(α) decreases with α, as expected. A linear dependence between the two was observed on a log-log scale (Fig. 1). According to the Cramer-Rao Lower Bound, the minimum variance of the α estimates is therefore expected to increase as the real parameter value increases. The real value of k has a more limited influence on IF(α) at the same scale (Fig. 1).

The variance of the α estimates obtained from simulated i.i.d. samples indeed increases with α (Fig. 2). Additionally, a positive bias is observed in the area of the parameter space delimited by low α values and high k values (i.e., α<30 and k>0.5 for N=1000), This situation is unlikely to occur in practice, as the strength of the coherent signal is expected to correlate with the number of scatterers per resolution cell (Fig. 3). In the rest of the parameter space, the bias is negative and becomes stronger as α increases and k decreases. This range of values is interpreted as a high number of scatterers per resolution cell, accompanied by the presence of some coherent signal. The negative bias is maintained even for k=0 (i.e. the absence of a coherent signal). According to [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF]Fig. 7], positive estimation biases are expected for k as α increases, therefore the underestimation of α appears to be at least in part correlated to the overestimation of k. One possible interpretation is the reduction in the number of effective scatterers resulting from the increasing regularity of their spatial arrangement and subsequent generation of a coherent signal [START_REF] Dutt | Ultrasound echo envelope analysis using a Homodyned K distribution signal model[END_REF]. Both the bias and variance decrease as the sample size increases.

B. Simulated Ultrasound Signals

Examples of envelope images corresponding to four concentrated simulated volumes, together with their histograms, are illustrated in Fig. 4. The fine differences in amplitude levels and their distributions are visible on the histograms, in the shift of the modes and heaviness of the tails.

Estimations conducted on data grouped as described in the previous section showed a monotonic increase of α with increasing scatterer numbers, but only for numbers lower than 43 per resolution cell (Fig. 5, Table 2), corresponding to the volume fraction 0.125. Beyond this limit, the estimates were no longer correlated to the number of scatterers per resolution cell. By comparison, the Nakagami parameter m reached a plateau between 10 and 20 scatterers/resolution cell, corresponding to volume fractions between 0.025 and 0.05. The variance reduction was visible as the number of samples increased, including for the largest sample size (50 000) where only 10 estimates were averaged. Additionally, the estimated mean values of α increased with increasing sample size, supporting the assumption made in the previous section that the negative estimation bias can be reduced by increasing the number of samples. In parallel to the statistical QUS study, the BSCs and the structure functions were also computed from the CREANUIS simulations, as detailed in the Appendix. A good agreement was found between BSCs and structure functions of simulated and theoretical data (see Fig. 8), which allows to demonstrate the ability of the simulations to reproduce the interference effects caused by spatial correlations among scatterer positions in concentrated media. Table 2. Means and standard deviations of the α estimates obtained from simulated signals (Fig. 5), for the interval where the mean estimates increase with the number of scatterers per resolution cell. N denotes the number of samples.

C. Tissue-mimicking Phantoms

The evolution of estimates obtained from phantoms was similar to the one observed in the case of the simulated signals. A monotonic increase of α with scatterer concentration was observed up to 40 scatterers/resolution cell (Fig. 6, Table 3), corresponding to a volume fraction of 0.1. Beyond this value, estimates and scatterer number were no longer correlated. The largest sample size (50 000) allowed for the best distinction between concentrations, although averaging was performed on only 4 estimations. As in the previous cases, the mean estimates of the α parameter achieve a better distinction of the concentrations than the mean estimates of the m parameter. However, because of the high variance, it would be more practical to represent a derived parameter with a lower variance. The small sample size (N=1000 samples) bias remained consistent with the one observed in the previous two simulations. No significant difference was observed between the larger sample sizes (N=10000 samples and N=50000 samples) in the experiment, presumably due to the small number of estimates that were averaged. 3. Means and standard deviations of the α estimates obtained from signals acquired from particle phantoms (Fig. 6), for the interval where the mean estimates increase with the number of scatterers per resolution cell. N denotes the number of samples.

IV. DISCUSSION

Three types of data were used to analyze the predictive value of the Homodyned K scatterer clustering parameter α for media concentration: simulated i.i.d. samples from distributions with imposed parameters, simulated ultrasound signals and ultrasound signals acquired experimentally from tissue-mimicking particle phantoms. Estimations performed on the simulated distributions are affected by high bias and variance, both of which increase with the true value of α, as implied by the decreasing Fisher Information. In [START_REF] Abraham | Novel physical interpretations of K-distributed reverberation[END_REF], the Fisher Information of the related K-distribution was computed and found to be inversely proportional to 1/α 2. for large values of α. A linear dependence between IF(α) computed in the present study and 1/α 2 was also observed, although slope values vary with k. IF(α) decreases with α linearly on a log-log scale (Fig. 1).

For ultrasound speckle, increasing α values are expected to correlate with increasing numbers of scatterers per resolution cell. Indeed, estimations performed on both simulated and experimental ultrasound data support this claim, up to a limit of approximately 40 scatterers per resolution cell. As for any statistical estimation, the mean and variance of the estimates are also affected by the sample size. The variances of the estimates from simulated i.i.d. samples (Fig. 2) agree with those of estimates from ultrasound signals, for the respective sample sizes and estimated mean values (Table 2 andTable 3). Moreover, a negative bias is expected when dealing with concentrated media, respectively with large numbers of scatterers per resolution cell. As the sample number increases, the absolute value of this negative bias is reduced, thus implying that the estimated values increase and approach the true value. The results obtained on ultrasound data are in agreement with this assumption. However, if we refer to section III-A for the expected biases when using the available number of samples, it can be considered that the actual values of α for the respective densities were not accurately obtained.

The parameter m of the Nakagami distribution was estimated throughout the study, for comparison purposes. For the concentrated media where the mean estimates of m saturated around the value of 1, the mean estimates of α achieved a better distinction between concentrations. Nevertheless, the estimation of α is increasingly unreliable as the number of scatterers increases and the distribution approaches a Rician [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF]. We found that the values of α relate consistently to the number of scatterers per resolution cell up to a limit of approximately 40 scatterers per resolution cell, corresponding to a volume fraction of 0.1 at a central frequency of 22 MHz. Both in simulation and in the experiment, the estimates of the two parameters (α and m) decrease again after reaching their peak values, Both in simulation and in the experiment, the estimates of the two parameters (α and m) decrease again after reaching their peak values, suggesting either model inadequacy or numerical problems beyond the respective concentrations.

The evolutions of the two studied parameters are similar to the evolution of BSCs measured in previous experiments conducted on dense suspensions of polyamide microspheres [START_REF] Franceschini | Experimental assessment of four ultrasound scattering models for characterizing concentrated tissue-mimicking phantoms[END_REF]. These experiments generally showed a good agreement between experimental BSCs and theoretical BSCs computed using the concentrated structure factor model in the 6-22 MHz frequency bandwidth [START_REF] Franceschini | Experimental assessment of four ultrasound scattering models for characterizing concentrated tissue-mimicking phantoms[END_REF]Figs. 6,7 and 8]). However, the structure factor model slightly overestimated the experimental BSC magnitude for volume fractions greater than 0.1 ([13, Fig. 8]). One may wonder whether the multiple scattering in our experiments might explain the discrepancy between theory and experiment. In [START_REF] Franceschini | Experimental assessment of four ultrasound scattering models for characterizing concentrated tissue-mimicking phantoms[END_REF], the length of the rectangular window used to analyze the data was equal to 1.3 mm, corresponding to 15 wavelengths at a center frequency of 17.5 MHz. In the current study, in order to mitigate potential effects of multiple scattering, the analysis was performed using a much shorter rectangular window of 0.3 mm, corresponding to approximately 4.5 wavelengths at a center frequency of 22 MHz. By considering the size and acoustic properties of polyamide particles for the highest volume fraction (25%), the Keller model [START_REF] Derode | Influence of correlations between scatterers on the attenuation of the coherent wave in a random medium[END_REF] [49] predicts a scattering mean free path of 0.95 mm at 22 MHz. Therefore, because the return distance of the ultrasound echo (0.6 mm) is inferior to the scattering mean free path, multiple scattering is negligible and does not explain the discrepancy between theory and experiments.

The experimental results were obtained after compensating the signals for attenuation (eq. 18) using the mean estimated attenuation coefficients αatt (Table 1). Compensating with the lower bound values of αatt led to an underestimation of α, while compensating with the upper bound values accomplished the opposite. The mean values of α changed by an average of ± 1%, ± 2.4% and ± 3.4% for the three samples sizes, respectively, while the variances were unchanged. Overall, the variance of the attenuation coefficients did not determine a statistically significant change in the α estimates.

We note that the configurations of the simulation and the experiment are slightly different due to the empirical adaptation of the resolution of the simulated probe to that of the experimental probe. It was practical to approximate the axial resolution of the single-element probe by slightly shifting the central frequency and assuming a symmetrical frequency bandwith.

The interpretation of α is indeed difficult because of the large estimation variance. However, a representation of 1/α versus scatterer numbers on a log-log scale (Fig. 7) shows a decreased variance and better overall separation of the different concentrations. In addition, the dependency of 1/α on the number of scatterers per resolution cell is approximately linear in this domain (similar to the results obtained in [START_REF] Dutt | Ultrasound echo envelope analysis using a Homodyned K distribution signal model[END_REF]). The estimates obtained from experiments appear to systematically have higher values than those obtained from simulations, possibly due to differences in the approximation methods of the three-dimensional resolution cell size (the resolutions of the simulated images was estimated by measurement, whereas the resolutions of the acquisition probe were taken from the datasheet).

Using large amounts of independent samples can be challenging due to the requirement that they belong to the same distribution. In practice, this can be achieved if the scanned medium is homogenous in terms of scatterer concentration and therefore expected to produce speckle with the same statistical properties at any location. The values estimated in section III-A can serve as guidelines when working with different sample sizes (i.e. expected bias and variance). For ultrasound data, averaging α estimates from sets of 10 000 i.i.d. samples delivered values that differentiate between the different studied scatterer numbers (up to 40 per resolution cell) with standard deviations of approximately 50% of the estimated means. Moreover, averaging only 10 estimates obtained on 10 000 i.i.d. samples (equivalent to a ROI of 1 cm in width and 0.5 mm in depth using the configuration of our 30 MHz simulation) yields approximately the same values as those presented in Tables 2 and3. It is important to notice that the simulation example, as well as the proposed experiment, are limited to small rectangular windows in depth (0.5 mm and 0.3 mm, respectively), corresponding to the length of the focal zone. In clinical practice, the use of a conventional multi-element probe would allow to deal with larger rectangular windows (thus increasing the number of i.i.d). Another way to increase the number of uncorrelated samples is to perform angular compounding as proposed by Hruska and Oelze [START_REF] Hruska | Improved parameter estimates based on the Homodyned K distribution[END_REF].

V. CONCLUSION

The present paper summarizes the set-up and results of an analysis of the scatterer concentration characterization power of the Homodyned K scatterer clustering parameter α. By using the parameter estimation method presented in [START_REF] Destrempes | Estimation method of the Homodyned K distribution based on the mean intensity and two log-moments[END_REF], a series of estimates were computed on simulated i.i.d. samples from Homodyned-K distributions with varying parameters, simulated ultrasound signals, and experimental ultrasound signals obtained from sparse and concentrated particle phantoms. The parameter α relates to the number of scatterers per resolution cell and can differentiate between scatterer concentrations especially in relatively sparse media. We find that its use is not limited to partially developed speckle, but can also extend to Rayleigh and Rician statistics if a large number of uncorrelated samples are available. We have successfully modeled scatterer volume fractions of up to 0.1 imaged at a frequency of 22 MHz by using estimates of α. We therefore suggest that it is possible to use α for the characterization of concentrated media such as blood or tumors with densely packed cells, provided that the acquisition of a sufficient amount of independent data samples is possible under similar experimental conditions.

APPENDIX

The appendix aims at comparing the BSCs and structure functions predicted by a model assuming identical non-overlapping spheres to the simulated BSCs and structure functions computed from CREANUIS simulations (presented in section II.C.). For each studied volume fraction, the simulated echoes were selected with a rectangular window of length 1 mm and the power spectra of the backscatter echoes were averaged to provide P. The simulated backscatter coefficient BSCsim was obtained using a reference scattering medium as follows ( [START_REF] Wang | An approach for measuring ultrasonic backscattering from biological tissues with focused transducers[END_REF]): BSC 567 ( ) = 8(9) 8 : (9) BSC ; ( ) [START_REF] Shankar | Classification of ultrasonic B-mode images of breast masses using Nakagami distribution[END_REF] where k is the wavenumber, P is the backscattered power spectrum from the tissue to be characterized, Pl is the backscattered power spectrum from the simulated volume having the lowest volume fraction (φl=0.005) and BSCl is the theoretical BSC for an ensemble of identical fluid spheres computed as BSC ; ( ) = < = 9 " > ? @ A B @ CD @ E F 9G H I (2 0)J [START_REF] Trop | The added value of statistical modeling of backscatter properties in the management of breast lesions at US[END_REF] where a is the scatterer radius (equal to 6.5 µm), Vs is the scatterer volume, nl is the scatterer number density equal to (nl=φl/Vs), γz is the relative impedance contrast between the scatterer and the surrounding medium (equal to 0.1) and j1 is the spherical Bessel function of the first kind of order 1. The simulated structure function was deduced from the two simulations with two different volume fractions: a lower volume fraction for a reference purpose and another volume fraction under consideration. The lower volume fraction (=0.005) is sufficiently low such that the structure function is assumed to be unity (incoherent scattering). The simulated structure function was obtained by [START_REF] Derode | Influence of correlations between scatterers on the attenuation of the coherent wave in a random medium[END_REF]:

567 ( ) = : KLM NOP (9) KLM Q (9) = : 8(9) 8 Q (9) (21

)
where n is the scatterer number density for the simulated medium under consideration. Fig. 8 shows the simulated BSCs and structure functions obtained from different volume fractions, as well as theoretical BSCs computed using the Structure Factor Model (SFM) [START_REF] Franceschini | Quantitative characterization of tissue microstructure in concentrated cell pellet biophantoms based on the structure factor model[END_REF]Fig. 8] and the theoretical structure functions based on [START_REF] Wertheim | Exact solution of the Percus-Yevick integral equation for hard spheres[END_REF], in solid lines. One can notice that as the volume fraction increases, the simulated structure function curves become smoother, as observed previously in [START_REF] Franceschini | Quantitative characterization of tissue microstructure in concentrated cell pellet biophantoms based on the structure factor model[END_REF]. A good agreement is present between simulated and theoretical BSCs and structure functions.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Fisher Information for increasing values of the scatterer clustering parameter α (obtained by numerical simulation). The values on the xaxis correspond to α values between 1 and 80, and the different curves correspond to different values of the coherent-to-diffuse signal ratio k.

Fig. 3 .

 3 Fig. 3. Biases for the parameter α as a function of the true values of α and k, for three sample sizes (N): 1 000, 10 000, 50 000. Each value was obtained by averaging 100 estimations.

Fig. 4 .

 4 Fig.4. (a) Extracts of envelope images obtained from simulated signals for 4 scatterer volume fractions: 0.1, 0.125, 0.15 and 0.2. (b) Envelope histograms corresponding to the same volume fractions.

Fig. 5 .

 5 Fig. 5. Results of the estimations performed on simulated ultrasound signals: (a) the Homodyned K parameter α and (b) Nakagami parameter m.The solid lines mark the interval where the mean estimates increase with the number of scatterers per resolution cell.

Fig. 6 .

 6 Fig. 6. Results of the estimations performed on signals acquired from particle phantoms: (a) the Homodyned K parameter α and (b) the Nakagami parameter m. The solid lines mark the interval where the mean estimates increase with the number of scatterers per resolution cell.

Fig. 7

 7 Fig.7 Log-log scale representation of estimates of the inverse α parameter from both simulated and experimental signals. Scatterers/resolution cell values correspond to the 1.7 -43 range in the linear domain.

Fig. 8 .

 8 Fig. 8. (a) Comparison between the simulated BSCs and the SFM predictions for several volume fractions. (b) Comparison between the simulated structure functions and the structure function of identical non-overlapping spheres. The solid lines represent theoretical values.

Table 1 .

 1 Equivalent number of scatterers per resolution cell (rounded) for the simulated phantoms. Equivalent number of scatterers per resolution cell (rounded) and measured total attenuation coefficients (αatt) for each of the 12 particle phantoms.

	Particle phantoms

ACKNOWLEGMENTS

This work was performed within the framework of the LABEX CELYA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). The authors would also like to thank Olivier Lombard for his useful inputs on result interpretation, as well as the anonymous reviewers who have provided an abundance of valuable questions and suggestions.