{'fraction_non_alphanumeric': 0.030588525703001256, 'fraction_numerical': 0.014214304815584309, 'mean_word_length': 5.08664079128584, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Multiparametric Volumetric Quantitative Ultrasound Imaging Technique for Soft Tissue Characterization

Quantitative ultrasound (QUS) offers a non-invasive and objective way to quantify tissue health. We recently presented a spatially adaptive regularization method for reconstruction of a single QUS parameter, limited to a two dimensional region. That proof-of-concept study showed that regularization using homogeneity prior improves the fundamental precision-resolution trade-off in QUS estimation. Based on the weighted regularization scheme, we now present a multiparametric 3D weighted QUS (3D QUS)imaging system, involving the reconstruction of three QUS parameters: attenuation coefficient estimate (ACE), integrated backscatter coefficient (IBC) and effective scatterer diameter (ESD). With the phantom studies, we demonstrate that our proposed method accurately reconstructs QUS parameters, resulting in high reconstruction contrast and therefore improved diagnostic utility. Additionally, the proposed method offers the ability to analyze the spatial distribution of QUS parameters in 3D, which allows for superior tissue characterization. We apply a three-dimensional total variation regularization method for the volumetric QUS reconstruction. The 3D regularization involving N planes results in a high QUS estimation precision, with an improvement of standard deviation over the theoretical rate achievable by compounding N independent realizations. In the in vivo liver study, we demonstrate the advantage of adopting a multiparametric approach over the single parametric counterpart, where a simple quadratic discriminant classifier using feature combination of three QUS parameters was able to attain a perfect classification performance to distinguish between normal and fatty liver cases.

Abstract. Quantitative ultrasound (QUS) offers a non-invasive and objective way to quantify tissue health. We recently presented a spatially adaptive regularization method for reconstruction of a single QUS parameter, limited to a two dimensional region. That proof-of-concept study showed that regularization using homogeneity prior improves the fundamental precision-resolution trade-off in QUS estimation. Based on the weighted regularization scheme, we now present a multiparametric 3D weighted QUS (3D QUS) imaging system, involving the reconstruction of three QUS parameters: attenuation coefficient estimate (ACE), integrated backscatter coefficient (IBC) and effective scatterer diameter (ESD). With the phantom studies, we demonstrate that our proposed method accurately reconstructs QUS parameters, resulting in high reconstruction contrast and therefore improved diagnostic utility. Additionally, the proposed method offers the ability to analyze the spatial distribution of QUS parameters in 3D, which allows for superior tissue characterization. We apply a three-dimensional total variation regularization method for the volumetric QUS reconstruction. The 3D regularization involving N planes results in a high QUS estimation precision, with an improvement of standard deviation over the theoretical 1/ √ N rate achievable by compounding N independent realizations. In the in vivo liver study, we demonstrate the advantage of adopting a multiparametric approach over the single parametric counterpart, where a simple quadratic discriminant classifier using feature combination of three QUS parameters was able to attain a perfect classification performance to distinguish between normal and fatty liver cases.

Introduction
The role of ultrasound, as a portable triage and monitoring tool, is becoming increasingly important, especially in the current context of a global pandemic. Conventional ultrasound images are formed from backscattered radio frequency echo signals created as a result of interaction with macro-and micro scale tissue structures. The B-mode image processing, including envelope detection and log compression, highlights the large-scale (>ultrasound wavelength) features such as the organ boundaries, while discarding the frequency-dependent smaller scale (<ultrasound wavelength) structural information. The speckle pattern and echogenecity of the B-mode image, though related to the complex organization of smaller scale tissue structure and their averaging scattering strength, render the US imaging qualitative, system-dependent and open to user interpretation. Three decades of research in the field of quantitative ultrasound (QUS), followed by the pioneering works by Insana et al. [20] and Lizzi et al. [8], established system and user independent techniques to quantify acoustic and microstructural properties by parametrizing the frequency dependent contents in the backscattered echo. Due to the reduced dependence on user interpretation as well as the efficacy for tissue typing and monitoring disease progression substantiated by extensive experimental research [29,36], QUS imaging promises to further enhance the clinical utility of ultrasound.
Spectral based parametric imaging forms a major branch of QUS, which includes measurement of attenuation coefficient estimate (ACE), backscatter coefficient and scatterer properties such as effective scatterer diameter (ESD) and effective acoustic concentration (EAC) of scatterers. ACE is a frequency dependent acoustic property that measures the ultrasound energy loss with propagation depth due to the combined effect of scattering and absorption. Backscatter coefficient, on the hand, measures the ultrasound intensity that is scattered in the backward direction (normalized ultrasound backscattered power per unit volume of scatterers). The integrated backscatter coefficient (IBC) is a measure of the backscatter strength, which is estimated by integrating the backscatter coefficient in the effective frequency band of the transducer. Backscatter coefficient can further be used to estimate the ESD of the dominant sub-resolution (<wavelength) scatterer, with the effective scale being determined by the incident ultrasound frequency. The spectral based QUS parameters have been successfully used in numerous clinical applications, including cancer detection in prostrate and lymph nodes [11], bone health assessment [24], cervical ripening detection [12], breast lesion characterization [32], and placenta characterization [9].
Liver steatosis quantification using QUS has undergone a recent surge in an effort to address the alarming increase in non-alcoholic fatty liver disease prevalence rate (25% globally) [28]. Previous works demonstrated that the cellular ballooning due to fat infiltration in the fatty liver affects the ultrasonic scattering process, resulting in an increase in ACE and Backscatter coefficient [4,40]. Several clinical studies on QUS showed promising results to quantify the fat content in liver, with efficacy comparable to magnetic resonance imaging pro-ton density fat fraction (MRI-PDFF) and biopsy, the current gold standard for steatosis grading [2,10,13,26].
Current challenges towards the clinical translation of QUS imaging include large estimation bias and variance due to tissue inhomogeneity and precisionresolution trade-off [27,37]. The resolution-precision trade-off is particularly critical for thin (e.g. dermis) and heterogeneous (e.g. placenta) tissue characterization. Available smaller windows result in noisy power spectral estimates due to limited spectral resolution and spatial variation noise inherent in ultrasound scattering. Deviation from the homogeneity assumption arising from the variation in underlying scatterer statistics [39], on the other hand, is likely to affect most biological tissue to different extent. Rubert et al. demonstrated that the high variance (with standard deviation being around 50% of the mean value) in ACE for ex vivo bovine liver is attributed to the deviation from the Rayleigh distribution, rather to the biological variability [44]. A large variance due to spatial heterogeneity was also reported for a cervix remodelling assessment study, where the difference in QUS parameters between early and late pregnant cervix were obscured by the intra-class variation [17,31]. The reliability of Controlled Attenuation Parameter (CAP, Fibroscan, Echosens, France), a promising commercial tool to detect and quantify hepatic steatosis, has been reported to be compromised due to the variability in attenuation estimates. For example, an IQR>40dBm was found to be indicative of poor diagnostic accuracy [52].
Large axial window lengths and averaging over lateral and elevational RF lines reduce the spatial variation noise due to random scatterer positioning. These approaches offer improved accuracy and precision, assuming uniform scattering within interrogated volume. Averaging also improves the power spectral estimation [27]. Studies show that the Welch method, a refined periodogram that estimates the averaged periodogram of overlapped subwindows, was superior compared to rectangular, Hamming and Hanning windowing [27]. While these techniques attain an improved estimation at the expense of spatial resolution, more recent works explore the possibility of extending the precisionresolution trade-off. Regularization incorporating a spatial prior was successful in simultaneously improving the resolution and the precision of QUS estimation for homogeneous regions [6,41,48]. However, uniform regularization, which was used in [6,41,48], might not be appropriate for biological tissue with the presence of inhomogeneities, as it will lead to oversmoothing in homogeneous regions in an attempt to compensate for the local inhomogeneities. In our previous work, we presented a proof-of-concept study based on a spatially weighted total variation regularization method, where the adaptive regularization parameter was a function of local inhomogeneity [10]. The study was limited to the measurement of a single QUS parameter, ACE, in a 2D region-of-interest and was validated using a small in-vivo liver study.
Another factor that might affect the clinical uptake of QUS imaging is the provision of visual guidance. CAP, which is performed in A-mode, lacks visual guidance, limiting ability to select the optimum region-of-interest. Several recent works on liver imaging have addressed the importance of B-mode and parametric image guidance to select an appropriate region-of-interest (ROI) that excludes vessels and extrahepatic areas [13,23]. 3D QUS imaging can further be advantageous over the 1-D or 2-D counterpart, allowing volumetric reconstruction, multimodal registration along with enhanced visualization and interpretation.
In this work, we present a multiparametric 3D weighted QUS (3D QUS) imaging system using a 3D total variation regularization method incorporating spatially variable regularization parameters. The proposed 3D adaptive total variation (TV) regularization exploits the spatial consistency information while preserving the variation in scattering properties in all three directions, thus would provide a rich spatial description of the tissue properties. The adaptive regularization is particularly important for QUS reconstruction in tissue regions with backscatter variation. We perform phantom validation, using phantoms with uniform and variable backscatter properties and demonstrate the superior QUS reconstruction performance of the proposed system compared to the baseline methods. Furthermore, the proposed method simultaneously generates multiple QUS maps, including ACE, IBC, and ESD, and therefore provides complementary information regarding the underlying tissue structure. The multiparametric approach could improve the class separability between normal and pathological tissue in the feature space. Toward this goal, we evaluate the benefit of a multiparametric approach over the single QUS counterpart in an in vivo liver study including healthy volunteers and patients with hepatic steatosis.

3D QUS: Proposed Volumetric Quantitative Ultrasound Parameter Estimation Algorithm
In our previous work [10], we presented a single QUS parameter estimation method limited to a 2D region-of-interest. In the problem formulation, we utilized a simplified model for the backscatter coefficient, ignoring the frequency dependent term. Herein, we introduce a 3D multiparametric QUS system, using an adaptive three-dimensional total variation regularization method (section 2.1). As the adaptive regularization parameter, we utilize a 3D spatially weighted matrix as a function of envelope SNR (signal-to-noise ratio) deviation (section 2.1). We further adopt an improved model of the backscatter coefficient considering its frequency dependence. Using this model, we formulate equations to extract scatterer statistics, such as effective scatterer diameter and effective acoustic concentration (section 2.1).
Attenuation Coefficient Estimate and Integrated Backscatter Coefficient Measurement A pulse-echo ultrasound system, with the assumption of weak scattering (i.e. the Born approximation) and negligible effect of diffraction and refraction, can be modeled as a simple frequency domain multiplication of pulse characteristics and backscatter characteristics [47]. The backscattered signal intensity, S(f, z) obtained from a time-gated radio frequency signal window centered at depth z from the transducer surface at frequency f can be written as: where G(f ) represents the transducer transfer function and transmit pulse. A(f, z) represents the total attenuation effect from the transducer surface to the center of the respective time-gated window at depth z and B(f ) is the backscatter term. For most soft tissue, A(f, z) is approximately linearly proportional to the frequency [49]: where α (Np/cm/MHz) is the effective amplitude attenuation coefficient estimate for the round trip of the depth z.
The backscatter term, B(f ), in its simplest form, can be expressed as a scalar parameter, as has been adopted in the scattering model used by Field II [22] and previous ACE computation works [7,10]. However, a more appropriate model is required to capture the non-linear frequency dependence of the BSC term [47,48]: The reference phantom method [53], a widely adopted method for QUS estimation, removes the system dependence (e.g. cancels out the system dependent term, G(f )) by normalizing the signal intensity backscattered from the tissue sample utilizing the one obtained from a well-characterized reference phantom under the same transducer and system settings. Therefore, the ratio of the intensity backscattered from the tissue sample to that from the reference at (i, j, k) location at frequency f p can be written as: where s and r denote sample under experiment and the reference, respectively. Taking the natural logarithm of equation 4 leads to Substituting the following variables in equation 5 as we obtain, y = φα + β + n ln f The system is defined in a N 1 × N 2 × N 3 (N = N 1 N 2 N 3 ) spatial grid and a frequency band discretized at M points is considered. We can represent the volumetric variables in vector format, for any variable a ∈ R N with elements , reducing the forward model in equation 7 to where Here, φ m ∈ R N ×N is a diagonal matrix, with the diagonal elements evaluated at propagation depth z ∈ R N for frequency f m (m ∈ [1, M ]) is given by The vector y ∈ R N M concatenates M sub-vectors y m ∈ R N , with each y m representing the volumetric power spectra ratio term evaluated at frequency f m and I ∈ R N ×N is a identity matrix. N (0, σ N ) is the Gaussian noise present in the measurement of y. QUS parameters α, β, and n can abruptly change at the interface between two different tissue types, but these parameters will undergo gradual change (i.e. sparse spatial variation) within each tissue type. Therefore, with the piecewise continuous assumption in all three spatial direction, we propose a 3D total variation regularization approach for the reconstruction of x = [α, β, n] from the noisy estimation y: where the first term is the data fidelity term, the last three terms are the anisotropic 3D total variation regularization terms, and λ 1 , λ 2 and λ 3 are the regularization weights. The 3D spatially weighted TV regularizer for a variable η is defined as where the 3D weighted matrix W i,j,k is uniform for α, and spatially variable for β and n. h i+1,i , h j+1,j , h k+1,k are distances between the adjacent axial, lateral and elevational pixels respectively at corresponding spatial positions. The TV regularizer [45] minimizes the sum of the gradients by imposing sparsity promoting l 1 -penalty in all three directions of the reconstructed QUS volume. By promoting signals with sparse gradients, TV minimization recovers piecewise-constant volumetric parameters. Solving (10) and (6) will give the values of volumetric α s , β s and n s . From the effective attenuation coefficient estimate, α obtained at depth z, the local attenuation coefficient estimate, α local can be computed as: The integrated backscatter coefficient (IBC) at (i, j, k) location can be computed as: Selection of the Spatially Weighted Regularization Parameter Spectralbased methods, such as the reference phantom method, make the assumption that RF signal generation is a stationary scattering process, arising from a large number of randomly distributed scatterers identical in size. The RF envelope of such speckle pattern (namely fully developed speckle) can be modelled using the Rayleigh distribution [15]. There will be large measurement error and variance in QUS estimation when the underlying assumption for Rayleigh distribution is violated. Such conditions occur for ROIs with varying backscatter properties or finite scatterer density (< 10 scatterers per resolution cell) [39,44]. Therefore an adaptive regularization parameter is required, which is sensitive to the conditions leading to the underlying assumption for Rayleigh distribution. Envelope SNR (SNR env ), defined as the ratio of the mean to the standard deviation of the RF signal envelope, is constant for a Rayleigh distribution and equals to SNR opt =1.91. The SNR env value deviates from the Rayleigh level of 1.91 for conditions involving finite scatterer density (< 10 scatterers per resolution cell) or reduced degree of randomness (due to the presence of periodicity or clustering). Therefore, SNR env can be used to identify media with low scatterer density, a source of large variance in ACE measurement [44]. Though SNR env cannot discriminate media with different backscatter characteristic for large scatterer densities (> 10 scatterers per resolution cell), it can identify the boundary of two such media, as the Rayleigh assumptions of identical diameter and random distribution do not hold at the boundary.
We define a normalized SNR env parameter, envelope SNR deviation (∆SNR env ) [9], given at a spatial location (i, j, k) by Previous works qualitatively demonstrated the association between ∆SNR env and backscatter variation from the visual inspection of inhomogeneity, lacking the quantitative description of BSC parameters [9,10]. In this work, we investigate the dependence of ∆SNR env on the quantitative backscatter parameters (i.e. IBC) and therefore validate the efficacy of ∆SNR env to detect the backscatter variation within the Region-of-interest (ROI).
We perform a feasibility analysis on RF data acquired from a custom built ultrasound phantom (Table 1) with three inclusions of similar ACE and different backscatter characteristic, manufactured by CIRS (Norfolk, VA, USA). Inclusion 1 (with low scatterer density) can be identified from the background with significantly high ∆SNR env , whereas inclusion 2 and inclusion 3 have similar ∆SNR env as the background (Fig.1b). However, the boundaries for all three inclusions are highlighted with increased ∆SNR env . This trend is also evident in the Rayleigh distribution fit to the RF envelope amplitude for selected ROIs, where the envelope amplitude in the ROIs on the boundaries deviate from the Rayleigh distribution (Fig. 1c). The spatial variation in ∆SNR env correlates with the change in decibel scaled IBC (dB) (Fig. 1d). Therefore, ∆SNR env is a suitable regularization parameter which can identify backscatter variation as well as sufficiently low scatterer density.
We propose a 3D spatially weighted matrix, {W i,j,k } ∈ R N as a function of ∆SNR env to adaptively regularize the backscatter parameters, n and β: where a and b are constants and ∆SNR env min is a nominal ∆SNR env value for homogeneous region. For ∆SNR env << ∆SNR env min , the weighting has little effect on the regularization. As the weight decreases with the increase in ∆SNR env , the regularization effect on the backscatter terms is relaxed.
Effective Scatterer Diameter (ESD) Measurement The backscatter coefficient, when expressed in terms of an acoustic intensity form factor, can be used to extract the effective scatterer diameter (ESD) and effective acoustic concentration (EAC). We utilize a Gaussian form factor F gauss = exp(−0.827k 2 a 2 ef f ) [21], which assumes scattering sources with continuously varying acoustic properties randomly suspended in a fluid-like media. The backscatter coefficient, B can be written as: Here, k = 2πf /c is the wave number, c is the sound speed, and a ef f is the effective scatterer radius. nγ 2 0 is the net scattering strength, which is the product of the concentration of scatterers n times the relative impedance difference between the scatterers and surrounding tissues,γ 2 0 . B(f ) (= βf n ) is the backscatter coefficient measured using the 3D adaptive regularization method (section 2.1).
Dividing by f 4 and then taking the natural logarithm on both sides of (14) we get: We can estimate a ef f and nγ 2 0 by fitting a regression line, The effective scatterer diameter (ESD = 2a ef f ) can be estimated from (16) and effective acoustic concentration (EAC = 10 log(nγ 2 0 )) can be obtained solving (16) and (17) [42]. The 'Model Inversion' technique [38] of measuring scattering properties by comparing measured backscatter spectrum to the theoretical one and fitting a regression line has been utilized in several previous works [34,42]. This technique is computationally efficient over the traditional minimum average squared deviation (MASD) technique [21]. However, both of   these techniques employ single time-gated window-based backscatter coefficient measures, which are prone to large measurement error. A previous work showed an improvement in QUS estimation using a nearest neighborhood averaging, which exploits the similarity among QUS measures of neighbouring windows [34].
In the proposed method, ESD is derived from a TV regularized backscatter measure and therefore exploits the sparse spatial variation of QUS parameters.

Data Sets
Phantoms  Table 1. For all the experiments, zone 1 of phantom 1 was used as the reference phantom.
In vivo Liver A cross-sectional study was performed, which included a cohort of twenty-one adult participants, with 13 healthy controls (defined as MRI-PDFF<5%) and 8 NAFLD patients (MRI-PDFF>5%). This study ( H14-01964) was performed under written informed consent after approval by the UBC Clin- . MRI-PDFF computed from the MRI data has been found to be strongly correlated with histological steatosis grading [35]. Therefore, MRI-PDFF was used as a gold standard for hepatic steatosis detection, where an MRI-PDFF> 5% was defined as presence of hepatic steatosis. The ultrasound data were acquired from the right lateral intercostal window of the volunteer, laid in a supine position on a shaker board propagating the mechanical excitation.

Data Acquisition and Processing
3D QUS System Overview The proposed volumetric QUS imaging technique was implemented on a standard 3D ultrasound system, that was also used for Shear Wave Absolute Vibro-Elastography (S-WAVE) developed in our lab [1]. The addition of steady state shear waves for S-WAVE has negligible effect (shown later) on the QUS measures so the data acquisition for QUS can be performed either simultaneously with S-WAVE or serially, i.e. after removing shear wave excitation. The 3D US system, used in this work, comprises an Ultrasonix Sonix-Touch ultrasound machine (Analogic, Richmond, BC, Canada), a m4DC7-3/40 4D transducer, and a linear voice coil actuator for shear wave motion generation. The transducer center transmit frequency was set to 3.33 MHz whereas the imaging depth was set to 15 cm with a focus at 13 cm. The customized motor control module allows the acquisition of volumetric RF data consisting of 15 planes spaced at 1 • angular distance, with 25 frames acquired at each plane location at an effective frame rate of 250 frame-per-second (fps). Each 2-D frame in the 3-D volume consists of 64 RF scanlines covering a lateral field-of-view of 79 • . The sampling rate was 20 MHz, resulting in 3870 samples in each scanline.
RF data processing and analysis The proposed method was implemented in MATLAB 2018a (The MathWorks Inc., Natick, MA, USA). The optimization problem was solved using the convex optimization toolbox CVX in MATLAB [16]. We compare the performance of the proposed 3D QUS method against the baseline reference phantom method (RPM) and a uniform 2D QUS method utilizing 2D uniform total variation regularization [6], applied on each of the elevational planes separately. We manually selected ROIs within the tissue or phantom, excluding the nearfield data suffering from ringdown and other artifacts associated with transducerto-phantom interface. Each ROI included 25 temporal frames. Each frame within the ROIs were divided into time-gated windows with 80% overlap in both lateral and axial direction. The optimum dimensions for the time-gated windows for RPM were found to be 20 scan lines laterally and 50 wavelength (λ) (600 samples or 2.3 cm) axially. Here, the axial extent of one wavelength is approximately 0.46 mm. On the other hand, the dimensions for the time-gated windows for the regularization methods were selected as 5 scanlines laterally and 1.2 cm (300 samples or 25λ) axially. The Welch method [51] was used to obtain the power spectrum from the RF scan lines within each time-gated window. We considered the −20 dB bandwidth of the received power spectrum as the usable frequency range.
The power spectra obtained from a single realization suffers from large variance due to statistical fluctuations arising from a random media. Compounding n statistically independent but equivalent realizations have been found to improve the power spectral estimate by reducing the standard deviation proportional to 1/ √ n. Previously, different spatial and angular compounding based on spatial translations or rotations have been adopted to decorrelate the signal allowing independent realizations of RF data [14,18,25]. In the proposed method, we perform averaging of the periodograms spatially and temporally, using 5 lateral scanlines and 25 temporal frames, to compute the power spectral density. To further improve the QUS estimation using RPM, we average the periodograms over the elevational planes resulting in a single QUS plane. The 2D QUS method applies 2D uniform regularization on each of the elevational planes separately, whereas the 3D QUS method applies 3D adaptive regularization using all the elevational planes.

Evaluation Metrics
The mean and standard deviation of ACE, IBC and ESD were measured and reported. The accuracy and precision were assessed using the mean absolute percentage error and standard deviation of absolute percentage error, respectively.
Absolute percentage error of a quantity a ∈ R N , is defined as: where AP E is the absolute percentage error, a estimated and a actual represent the reconstructed and the ground truth measure of a QUS quantity, respectively. The detectability of an inclusion in a QUS map was evaluated using contrastto-noise ratio (CNR) ratio: where µ and σ are the mean and standard deviation of a QUS quantity within the inclusion and the background denoted by the subscripts inc and bg, respectively. The inclusions were delineated from the B-mode images, refined using the position and dimension information provided by the manufacturer.
To evaluate steatosis detection performance for in vivo liver, the correlation between MRI-PDFF and QUS measures were calculated using Spearman's rank correlation coefficient. Also, the classification accuracy and the area under the ROC curve (AUROC) were measured to evaluate the classification performance using different combination of QUS parameters.

Phantom Validation
The ACE reconstruction results for phantom 2 (variable ACE and similar backscatter and scatterer size properties) are shown in Fig. 3. From the ∆SNR env map, it is evident that the variation in ACE does not cause deviation of SNR env values from SNR opt . With ∆SNR env < ∆SNR env min , the spatial weighting has little effect on the regularization. Therefore, QUS reconstruction results from the 2D QUS method with uniform regularization and the 3D QUS with adaptive regularization are similar, where both methods clearly identify the low attenuation inclusion from the high attenuation background. The RPM, on the other hand, detects the variation in ACE with a poorer resolution compared to the regularization methods. Figure 4 shows the results for phantom 3 (approximately uniform ACE and variable backscatter and scatterer size properties). As previously discussed in section 2.1, the variation in backscatter characteristic is associated with a variation in ∆SNR env map. The 3D QUS method, that incorporates the ∆SNR env information for enabling adaptive regularization, yields more accurate reconstruction of the the IBC map compared to its uniform regularization counterpart, 2D QUS ( Fig. 4(a)). The RPM failed to generate a reasonable IBC result. The ACE maps obtained from the three methods have also been shown in Fig. 4 (c). Both the RPM and the 2D QUS method exhibit consecutive undershoot and overshoot in ACE measurement at the boundary between medium with different backscatter properties, a trend described in previous literature [10,39]. The proposed 3D QUS method, on the contrary, was able to overcome the issue, yielding accurate ACE results.  The ESD map for phantom 3 has been shown in Fig. 5. The RPM and the 2D QUS method using uniform regularization were not successful to reconstruct the ESD map, as were expected from the inaccurate estimation of IBC map. The proposed 3D QUS method results in improved ESD reconstruction.
For the quantitative comparison of the QUS parameters, we compare the mean and standard deviation of the ACE in phantom 2 and IBC and ESD in phantom 3 for the inclusions and the background. For IBC, we compute the standard deviation of the IBC error expressed in dB. The results are shown in Fig. 6. We also report the mean absolute percentage error and standard deviation of absolute percentage error results in Table 2 and 3. In case of ACE, both the 2D QUS method and the 3D QUS method reach close to the ground truth value with high precision, enabling accurate ACE characterization. However, for IBC and ESD estimation, the proposed 3D QUS method consistently outperforms both the other methods.  We evaluate the CNR of the QUS parameters for different methods. Figure  7 shows the mean and standard deviation computed over the three inclusions in respective phantom. The mean CNR values obtained using the proposed 3D QUS method were larger than the other two methods. Compared to the baseline RPM, the 3D QUS method results in around 4.2 fold, 1.8 fold and 2.9 fold increase in the CNR for ACE, IBC, and ESD, respectively. Also, the 3D QUS improves the CNR by 1.4 times and 1.9 times compared to the 2D QUS method (uniform regularization) for IBC and ESD, respectively.

Effect of three-dimensional regularization
Unlike conventional QUS system, the proposed system in this work acquires multiple elevational planes of RF data. We perform a three-dimensional adaptive TV regularization to obtain the volumetric QUS maps. To study the effect  of three dimensional regularization, we vary the number of planes and evaluate the precision by computing standard deviation of absolute percentage error for ACE, IBC and ESD measures. We also investigate whether three dimensional regularization has any advantage over the two-dimensional counterpart. Therefore, we perform 2-dimensional adaptive regularization after averaging the estimated power spectra from multiple planes. The results have been shown in Fig. 8. Also, the 1/ √ N curves for IBC and ESD have been superimposed on the graphs, where N represents the number of planes used in the computation. For uncorrelated planes, the standard deviation of absolute percentage error curves should follow the 1/ √ N curve. The ACE variance remains unaffected (standard deviation of absolute percentage error¡ 0.5%) as the number of planes increases for this particular example of homogeneous phantom. For IBC and  ESD measures, increasing the number of planes improve the estimation variance for both 2D and 3D regularization methods. However, the 3D regularization not only outperforms the 2D regularization estimation precision, but the reduction in standard deviation of absolute percentage error measure improves over the 1/ √ N rate, which is the theoretical rate obtained by averaging N uncorrelated realizations.   In vivo Liver Fig. 9: Correlation coefficient between two 2-D radio-frequency (RF) signals separated by ∆N planes for a phantom and an in vivo liver example. Figure 9 shows that the radio-frequency data obtained from in vivo liver are more decorrelated compared to the one obtained from the phantom. The improved decorrelation between subsequent planes is an effective indicator of improvement in QUS measurement for in vivo liver compared to phantom.

Effect of Spatial Weighting in Total Variation Regularization
The phantom study demonstrates that the spatial weighting in the proposed 3D QUS method improves the QUS estimation, even in regions with high backscatter variation. In this section, we investigate the effect of spatial weighting on the selection of regularization parameter and QUS estimation performance. In 10(a) and (b), we compare the cost function value for different values of regularization parameters obtained at the ground truth QUS measures for the 3D adaptive QUS and 2D uniform QUS method. As it is difficult to visualize a 4-dimensional cost function, we set the values of λ 2 and λ 3 equal. We found that the spatial weighting keeps the optimum cost comparatively less variable for different λ 2 and λ 3 values for variable ACE region-of-interest (homogeneous as dictated by low ∆SNR env ). The effect is more pronounced for phantom 3 with variable backscatter property (high ∆SNR env ). Figure 10(c) demonstrates that the error in QUS parameter estimation can be optimized by choosing a suitable regularization parameter set for both the 3D QUS and the 2D QUS method. However, for the 3D QUS, the optimum ranges for regularization parameter coincide, whereas for the 2D QUS with uniform regularization, the ranges do not overlap. Therefore, it is possible to select an optimum regularization parameter for the 3D QUS method to optimize QUS parameter estimation irrespective of homogeneity (backscatter variation). On the contrary, the 2D QUS with uniform regularization would require different parameters set depending on the underlying tissue characteristic, which would not be feasible for clinical application.

In vivo Liver: Steatosis Detection Performance
The correlation results between MRI-PDFF and QUS measures have been displayed in 10. The ACE and IBC measures demonstrated significant positive associations with MRI-PDFF, whereas ESD was negatively associated to MRI-PDFF. The proposed 3D QUS method yields stronger correlation between ACE and MRI-PDFF (r = 0.877, p < .001) compared to the RPM (r = 0.415, p = 0.062) and the 2D QUS (r = 0.645, p = 0.002) method. The IBC measure had similar and significantly strong association with MRI-PDFF for both the 2D  From the correlation analysis, the 2D QUS and the 3D QUS method exhibit comparable performance. However, a feature analysis plot including all three QUS features in Fig. 11 shows that the 3D QUS results in an improved separability between the steatosis and non-steatosis group compared to the 2D QUS method.

2D QUS 3D QUS
We utilize a quadratic discriminant analysis (QDA) classification to classify the cases into steatosis and non-steatosis group. The classification accuracy and the area under the ROC curve (AUROC) were reported in Table 4. IBC, as a single feature, exhibits superior ability to discriminate while using the regularization methods: 2D QUS and 3D QUS. Among different feature combination, ACE and IBC, ACE and ESD, and all three feature combination yield comparable performance.

Discussion
Adopting a total variation regularization scheme extends the trade-off between estimation accuracy, precision, and spatial resolution. A window dimension of 30 scanlines (10 uncorrelated scanlines) and 35λ (13 pulse lengths) was found to be the minimum reliable dimension without incorporating regularization. This findings lie within the range of 7-15 pulse lengths and 10-20 uncorrelated scanlines, which are presented as the dimensions offering the best trade-off between spatial resolution, accuracy, and precision [43] for QUS estimation. In contrast, incorporating the total variation regularization into the 2D and the 3D QUS methods results in a dramatic improvement in precision of both ACE and IBC, where standard deviation of absolute percentage error remains below 5% for a window dimension of 10λ ( 3.75 pulse lengths) and 25λ ( 9.40 pulse lengths) axially and 5 ( 2 uncorrelated) scanlines, laterally. The mean of the absolute error of ACE and IBC remained below 6% and 40%, respectively for both the regularization methods. We can compare the accuracy performance with the previous work [48] using a regularization method, which was able to achieve a mean percentage error (MPE) of 0.46% in ACE and 327% MPE and 19.7% MPE in estimation of backscatter parameters, b and n. Our results remain consistent with the previous studies [50], where backscatter estimations have been found to suffer from significant measurement error, ranging within 2 orders of magnitude.
Both regularization methods yield similar performance for the homogeneous phantom (phantom 1) and variable ACE phantom (phantom 2), as can be seen in Fig. 3 and Fig. 4. With a similar scatterer properties (including IBC and ESD), the ∆SNR env values of both of these phantoms remain below ∆SNR env min , with little effect of spatial weighting. Therefore, the uniform QUS and the weighted QUS methods are essentially similar when the tissue under experiment is homo-geneous in terms of backscatter properties. The 3D QUS method outperforms the 2D QUS when the underlying tissue has variable backscatter properties, as is the case of phantom 3. Figure 4 shows that the ∆SNR env map was able to identify the boundaries where backscatter properties vary. Incorporating the ∆SNR env information in the 3D QUS enables adaptive regularization, and therefore results in more accurate reconstruction of the ACE and IBC map compared to the uniform regularization in 2D QUS. It should be noted that both the RPM and the 2D QUS method exhibit an underestimation and an overestimation of ACE centring the boundary between medium with different backscatter properties, a trend described in previous literature [10,39]. The 3D QUS method was able to overcome the issue, yielding accurate ACE and IBC results with close agreement to the values reported by the manufacturer.
For ESD computation, a Gaussian form factor has been used as it can successfully model the scattering properties of many soft tissues, including liver [20,33]. However, a spherical shell form factor is more appropriate for spherical glass beads, as were used in the phantoms in our experiments [19,20]. The Gaussian form factor would consistently overestimate the sphere size, resulting in an effective diameter greater than the physical diameter [19,27]. The effective diameter of the glass beads can be found by comparing the Gaussian form factor to the spherical shell form factor and computing the conversion factor for the particular physical diameter within the usable frequency range [27]. The physical diameters for inclusion 1, inclusion 2, inclusion 3 and the background medium in phantom 3 are: 220 µm, 120 µm, 40 µm, and 70 µm, respectively. With conversion factors 1.27, 1.31, 1.56, and 1.28, respectively, the effective diameters have been found to be : 342, 158, 51, and 90 µm. The accuracy of ESD estimation depends on the estimation accuracy of the backscatter parameter, n. Therefore, the accurate backscatter parameter estimation obtained by the 3D QUS further is translated to accurate ESD estimation as can be seen in Fig. 5.
The CIRS phantom was specifically designed to cover the wide range of QUS parameters to be expected in human liver in vivo with different stage of steatosis. Figure 6, it can be seen that the 2D QUS method and the 3D QUS method both reach close to the ground truth value with high precision for the entire range of ACE from 0.37 to 1.18 dB/cm/MHz, while both the mean and standard deviation of absolute percentage error remain below 10%. For IBC, the deviation of IBC from the ground truth values range from 0.90 to 3.63 dB using the 3D QUS method, where the mean absolute percentage error is below 25% and standard deviation of absolute percentage error is below 11%.
The accuracy and precision of ESD estimation further depends on the range of ka ef f , the product of the wave number k and the effective scatterer radius a ef f , which should be within 0.5 to 1.2 for optimum results [20]. For the frequency range of 2 MHz to 3.5 MHz used in our experiments, the ka ef f values for the scatterers in inclusion 1, 2, 3 and the background are: 1.40−2.44, 0.64−1.13, 0.21 − 0.36, and 0.37 − 0.64, respectively. Accordingly, we obtain high estimation error and variance for inclusion 1 and 3, as the corresponding ka ef f values are well outside of the optimal range. The high variance in the ESD estimation for background scatterers can partly be attributed to the boundary artifacts. Finally, the ka ef f for the scatterers in inclusion 2 is within the optimal range, resulting in an accurate and precise ESD estimation (Fig. 5, Fig.6c and Table  3).
The advantage of using spatial weighting in the proposed 3D QUS over the uniform QUS and the baseline RPM is further evident from the CNR study. The 3D QUS improves the contrast significantly, with an average CNR improvement of 316%, 79%, and 188% for ACE, IBC and ESD, respectively compared to the baseline RPM. Compared to the 2D QUS, the CNR improvement achieved using the 3D QUS method was 16%, 37% and 92% for ACE, IBC and ESD, respectively. This improved CNR would translate into better identification of abnormality in clinical application.
The proposed QUS imaging system utilizes the volumetric RF data acquired by a motorized 3D transducer. The use of a 3D system allows the volumetric QUS reconstruction, using a three-dimensional total variation regularization. Figure 8 shows that 3D regularization improves the estimation variance at a faster rate than the theoretical 1/ √ N curve obtained from averaging N uncorrelated realizations. The advantage of utilizing volumetric data could be further enhanced in case of in vivo liver data. The correlation between two radio-frequency signal data separated by ∆N planes, gives a measure of effectiveness of spatial compounding [18], where a low normalized correlation coefficient value (< 0.2) indicates independent data realizations. Figure 9 shows that the radio-frequency data obtained from in vivo liver get decorrelated at a faster rate with the increasing plane separation compared to the one obtained from the phantom. As the ground truth values for QUS variance in in vivo liver is unknown, we cannot obtain the exact measures of improvement with increasing data realizations. Nevertheless, the improved decorrelation between subsequent planes is an effective indicator of improvement in QUS measurement for in vivo liver compared to phantom.
The proposed 3D QUS method can be realized on a standard 3D ultrasound imaging system. This work was done as a part of a study together with S-WAVE. In a separate examination of data that repeated tests with and without external vibration for S-WAVE, the variability in QUS estimates were not found to be significantly different. Future work could further compare serial versus simultaneous S-WAVE and QUS acquisition.
Our QUS estimation results for liver in vivo were comparable to the MRI-PDFF performance, the current gold standard. The presence of hepatic steatosis was defined as MRI-PDFF> 5% [5]. From the correlation analysis, it was found that the ACE and IBC obtained using the proposed method have strong positive association with MRI-PDFF, whereas the ESD was found to be negatively correlated to MRI-PDFF (Fig. 12). The ACE and IBC correlation results are in agreement with previous studies on human liver in vivo [3,46], additionally with an improved correlation. On the other hand, this is the first work to report the correlation between ESD of human liver and MRI-PDFF. Nevertheless, previous studies on rabbit fatty liver found negative correlation between ESD and fat content of the liver [33]. Using different combination of these three QUS parameters estimated using the proposed method with a simple QDA classifier, we were able to discriminate between the steatosis and the non-steatosis class with an accuracy up to 100% and area under the receiver operating characteristic (AUROC) of 1. Among the QUS features, ESD was least effective to classify the two groups. However, the high contrast achievable using ESD (Fig. 7) suggests that ESD, though exhibit inferior performance to identify steatosis which is a diffuse condition, might be effective in characterizing focal diseases.
Using Youden index, the optimal threshold for ACE, IBC and ESD to detect the presence of steatosis were found to be 0.77 dB/cm/mHz, 8.5e −4 1/(sr.cm), and 146.3 µm. The ACE and IBC cut-off agree with the previous findings, where the ACE and IBC cut-off values were equal to 0.76 dB/cm/MHz and 2.7e −3 1/(sr.cm) [3]. Among the limitations of our study is the ESD values obtained for steatosis participants, which partially lie within non-optimal ka ef f < 0.5 range. Future work will involve application of comparatively higher frequency range to mitigate the problem. Also, we observe high variance in the ESD estimation for the liver study, which can be due to the presence of coherent component in the backscatter RF data. Therefore, we will further explore different signal decomposition techniques to separate the diffuse component from the coherent component [34]. Also, instead of using a Gaussian form factor, tissue-specific form factors can be derived from the impedance maps of underlying tissue and could be used for improved ESD computation [30]. Another limitation of the in vivo liver study is that the MRI-PDFF images and the ultrasound images are not registered. Though care has been taken to ensure that the planes-of-interest in MRI and ultrasound roughly coincide, a registration scheme would allow direct comparison with possible identification of focal lesions.

Conclusion
In this paper, we present a multiparametric 3D weighted QUS imaging system, which simultaneously measures three biomarkers for tissue characterization including ACE, IBC and ESD. The method incorporating a spatially adaptive regularization was able to achieve improved precision-resolution trade-off in QUS reconstruction in the presence of underlying tissue homogeneity. The threedimensional regularization further improves the QUS estimation precision with a standard deviation reduction faster than the theoretical 1/ √ N rate achievable by averaging of N uncorrelated realizations. With strong correlation to the current gold standard MRI-PDFF, each of the three parameters demonstrate the potential to detect hepatic steatosis. The combination of the QUS parameters further improve the steatosis classification performance. The proposed volumetric QUS system, built within an elastography system S-WAVE, therfore promises a comprehensive assessment tool for non-alcoholic fatty liver disease.