A deep learning approach for patchless estimation of ultrasound quantitative parametric image with uncertainty measurement

TL;DR

This paper introduces a deep learning method to estimate ultrasound quantitative parametric images, specifically scatterer number density, without the need for patch-based analysis, while also providing uncertainty measurements for the predictions.

Contribution

The study presents a novel deep learning approach that eliminates the patching process in QUS parameter estimation and incorporates uncertainty quantification.

Findings

Accurate estimation of scatterer number density without patching.

Uncertainty maps effectively indicate confidence levels in predictions.

Method outperforms traditional patch-based estimation techniques.

Abstract

Quantitative ultrasound (QUS) aims to find properties of scatterers which are related to the tissue microstructure. Among different QUS parameters, scatterer number density has been found to be a reliable biomarker for detecting different abnormalities. The homodyned K-distribution (HK-distribution) is a model for the probability density function of the ultrasound echo amplitude that can model different scattering scenarios but requires a large number of samples to be estimated reliably. Parametric images of HK-distribution parameters can be formed by dividing the envelope data into small overlapping patches and estimating parameters within the patches independently. This approach imposes two limiting constraints, the HK-distribution parameters are assumed to be constant within each patch, and each patch requires enough independent samples. In order to mitigate those problems, we employ a deep learning approach to estimate parametric images of scatterer number density (related to HK-distribution shape parameter) without patching. Furthermore, an uncertainty map of the network's prediction is quantified to provide insight about the confidence of the network about the estimated HK parameter values.

Tables

Table 1. Table 1: The simulation results of the proposed method.

	RMSE	RRMSE	MAE
Patchless CNN	1.12$\pm$0.59	0.42$\pm$0.19	0.85$\pm$0.55

Equations

$$P_{HK}(A|\epsilon ,{\sigma }^2,\alpha )=A{\int }_0^{\infty }u{J}_0(uϵ)J_0(uA)(1+\frac{u^2{\sigma }^2}{2}{)}^{-\alpha }du$$

$$s_{(a,l)}=TR{F}_{(a,l)}\ast h_{(a,l)}$$

$$h[a,l]=e^{-\frac{1}{2}(\frac{a^2}{{{\sigma }_a}^2}+\frac{l^2}{{{\sigma }_l}^2})}\times cos(2\pi f_{c}a)$$

$$\begin{array}{c}X=<Ilog(I)>/<I>-<log(I)>,\\ U=<log(I)>-log(<I>),\end{array}$$

$$\begin{array}{c}\tilde{S}=\frac{1}{N_f}\sum _{i=1}^{N_f}\widetilde{S_i},\\ Uncertainty=\frac{\sqrt{\frac{1}{N_f}{\sum }_{i=1}^{N_f}(\widetilde{S_i}-\tilde{S}{)}^2}}{\tilde{S}},\end{array}$$

$$loss=\frac{1}{N}\sum (lo{g}_{10}(\tilde{S})-lo{g}_{10}(S){)}^2$$

Full text

1 INTRODUCTION

Quantitative ultrasound (QUS) finds scatterers’ properties which are highly related to the tissue microstructure [1]. Different QUS parameters such as speed of sound, scatterer number density, and backscattering coefficient have been employed for characterization of tissues. Among different QUS parameters, scatterer number density has been found to be useful for different clinical applications such as liver fibrosis detection [2] and Hepatic steatosis assessment [3]. Scatterer number density can be quantified by modeling the ultrasound echo amplitude using different probability density functions. The homodyned K-distribution (HK-distribution) is a well-known distribution that has been widely used to quantify the scatterer number density and characterizing the tissues. This distribution can model diverse scattering scenarios but requires a large number of independent samples [4].

The HK-distribution parametric images are formed by dividing the envelope data into small overlapping patches and estimating parameters within the patches independently. Estimating the parametric images is challenging and the parametric images are often noisy but still have been found useful in clinical applications to detect different abnormalities [5]. The main challenge in estimating the HK-distribution parametric images is the small size of the patches. Increasing the size of the patches might not be helpful since the spatial resolution will be lost and the heterogeneity inside the patch is increased (due to heterogeneity of tissue types or spatially variant nature of point spread function). Correlated samples also should be skipped since they introduce bias to the estimation of HK-distribution parameters. Therefore, a large number of samples inside each patch should be ignored to reduce the correlation between the samples which further reduces the number of samples.

In this paper, an ultrasound simulation method is employed to generate the training data. Unlike sampling from HK-distribution, this method of simulation data generation contains correlated samples which is a more realistic method and closer to the real ultrasound data [6]. In addition to this, multiple frames from experimental phantom data are employed to obtain a more accurate parametric images and quantify the uncertainty map of the network. The uncertainty map would enable the clinicians to find out how reliable the estimations are in each region.

2 Materials and Methods

2.1 HK-distribution

The HK-distribution can be formulated as [4]:

[TABLE]

where $A$ denotes the envelope of the backscattered echo ultrasound data, $u$ is the variable of integral that needs to be integrated from zero to infinity, $J_{0}(.)$ is the zero-order Bessel function, and $\alpha$ is the scatterer clustering parameter which reflects to the scatterer number density. The parameter $\alpha$ and the ratio of coherent signal power ($\epsilon^{2}$) to the diffuse one ($2\sigma^{2}\alpha$) denoted as $k$ have been employed as the parameters of HK-distribution to characterize the tissues.

2.2 Data Generation

A diverse dataset with known scatterer number density is required to train the network. We followed the fast grid-based method [6] with the difference that here the scatterer number density can be any value in the range of 1-20. Assuming weak scattering, the ultrasound RF data can be obtained by the 2D convolution of Point Spread Function (PSF) and the Tissue Reflectivity Function (TRF) [7, 6].

[TABLE]

Where $h$ denotes the PSF. In order to have a known scatterer number density, the PSF is assumed to be spatially invariant which can be described as:

[TABLE]

where the PSF is modeled as a 2D Gaussian modulated with a cosine with the center frequency of $f_{c}$ in the axial direction. The parameters ${\sigma_{a}}^{2}$, and ${\sigma_{l}}^{2}$ denote the axial and lateral width of the Gaussian function and they are related to the resolution cell. The TRF is the 2D map that contains all point scatterers. In each grid point, only one scatterer is allowed to be present and its amplitude is sampled from normal distribution. The mean of the normal distribution is randomly selected from values 1 to 5, the variance is fixed to 0.02 of the mean value. The TRF is constructed using shapes having different values of scatterer number density and mean scattering amplitudes. One sample generated by this method is shown in Fig. 1. We generated 11000 samples to train the network. It should be noted that the ground truth value of scatterer number density is known and depends on how the resolution cell size is defined. Here, we define the resolution cell in each direction as $3\times\sigma_{x}$, where x is the axial or lateral directions. Different definition of the resolution cell results in different value of the scatterer number density. Please refer to [6] for more information about the grid-based simulation method. In order to reduce the correlation, samples are skipped from both axial and lateral directions to have an image of size $256\times 128$, and an average correlation between samples of 0.28 is obtained.

2.3 Conventional Patch-based Methods

Hruska et al. employed SNR, skewness, and kurtosis to estimate HK-distribution parameters [8]. Destrempes et al. proposed XU estimator in which two log moment statistics (named as X and U) were employed, and reported improved parameter estimation [4]. The X, and U can be defined as:

[TABLE]

where $<.>$ denotes sample mean, and $I=A^{2}$ is the intensity of backscattered signal. The XU estimator iteratively solves a constraint optimization using bisection method to find the parameters. We employed this method for comparison.

2.4 Experimental Phantom Data

Several experimental phantoms have been employed to evaluate the method.

1. Layered Phantom: A phantom with a middle layer having different properties than the top and bottom layers reported in previous publication [9] has been employed. The phantom was constructed by an emulsion of ultrafiltered milk and water-based gelatin, and 5–43 $\mu m$ diameter glass beads were utilized as the source of scattering. A Siemens Acuson S2000 scanner (Siemens Medical Solutions USA, Inc.) with 18L6 probe linear transducer having the center frequency of 8.9 MHz was used for data collection. The layer with higher intensity has a higher backscattering coefficient and scatterer concentration. The backscattering coefficient of this layer is $6.37\times 10^{-3}$ $cm^{-1}sr^{-1}$ and it is $3.52\times 10^{-3}$ $cm^{-1}$ $sr^{-1}$ for other parts at the center frequency.

2. CIRS Phantom: Data from a multipurpose CIRS phantom (model 040GSE, Norfolk, VA, USA) was collected at Concordia university using E-CUBE 12 Alpinion machine by an L3-12H transducer with a sampling rate of 40 MHz and center frequency of 8.5 MHz.

3. Gammex Phantom: Data was collected at university of Wisconsin using a Verasonics Vantage 128 System (Verasonics,Kirkland, WA, USA) machine with an L11-5v transducer operating at 8 MHz.

2.5 Pre-Processing of experimental phantom data

Experimental phantom data may contain variations of intensity due to beamforming, focusing and time gain control (TGC) which are not present in the simulation data and leads to performance decay of parameter estimation. Inspired by the reference phantom method which cancels out the system effect [10], we employ an intensity normalization technique to reduce the system effect on the envelope data intensity. First, the average of envelope amplitudes in each depth is computed across several frames of uniform reference phantoms imaged by the same imaging settings. In the next step, a curve is fitted to the obtained average envelope amplitude. The inverse of the fitted curve can be utilized to normalize the experimental phantom data imaged by the same settings. The procedure is depicted in the Fig. 2. It can be observed that the normalized data has a more uniform intensity across different depth compared to the pre-normalized one.

2.6 Multi-frame acquisition and data processing for experimental phantoms

For the layered phantom, the probe is moved in the out-of-plane direction to collect 12 frames from this phantom. We skipped 7 samples in the axial and 1 sample in the lateral direction from the acquired envelope data to reduce the correlation from 0.75 (no skipping) to 0.141. The patch size for XU algorithm is $4.5\times 4.5$ $mm$ with $75\%$ overlap.

For the CIRS phantom, the probe is moved in the out-of-plane direction to collect 10 frames from this phantom. The patch size for XU algorithm is $4.5\times 4.5$ $mm$ with $75\%$ overlap. We skipped 4 samples in the axial and 0 sample in the lateral direction from the acquired envelope data to reduce the correlation from 0.95 (no skipping) to 0.49. In another experiment, we skipped 5 samples in the axial and 1 sample in the lateral direction to have the correlation of 0.36 to see the effect of different correlations. It should be noted that skipping more samples results in reduction of the bias but lower number of samples would be available to estimate the scatterer number density which makes estimation of small parts challenging.

For the Gammex phantom, the probe is moved in the out-of-plane direction to collect 12 frames from this phantom. The patch size for XU algorithm is $4.0\times 4.0$ $mm$ with $75\%$ overlap. We skipped 6 samples in the axial and 1 sample in the lateral direction from the acquired envelope data to reduce the correlation from 0.95 (no skipping) to 0.55. In another experiment, we skipped 13 samples in the axial and 1 sample in the lateral direction to have the correlation of 0.19 to see the effect of different correlations.

2.7 Prediction and uncertainty quantification

For experimental phantom data, several frames are collected by sweeping the probe in the out-of-plane direction. The final estimate and its corresponding uncertainty can be quantified as:

[TABLE]

where the final estimate ($\widetilde{S}$) is calculated by averaging the estimates ($\widetilde{S_{i}}$) obtained from $N_{f}$ frames. The proposed uncertainty measures how different the estimates are across different frames. If the estimates from all frames are similar, the uncertainty will be low. On the contrary, high uncertainty is expected when the estimates are different. It should be noted that this is a frame-wise uncertainty measurement meaning that if there is no difference in the estimated value across different frames, the uncertainty is low even if the estimated values are incorrect.

2.8 Network architecture and training

We employed Deeplab V3 which has shown excellent performance in semantic segmentation [11]. The last layer of the model is altered to have one output as the scatterer number density. The following loss function is employed:

[TABLE]

where $\widetilde{S}$ and $S$ are the predicted and ground truth scatterer number density, respectively. The log compression is employed to avoid the bias toward high scatterer number density values. Rectified Adam (RAdam) is employed as the optimizer which has shown outperforming Adam and be robust to the selection of learning rate [12]. The network is trained on single NVIDIA RTX 3090 with 24 GB of memory.

3 Results

3.1 Simulation Results

The method is evaluated on 1000 test data generated in a similar way of training data. Due to the fact that the value of the parameter $\alpha$ is not known for the simulation data, we cannot compare it with the XU estimator. The results are given in Table 1. The metrics are evaluated on each sample of test dataset and, the average and standard deviation are reported. Four samples of the simulation results are depicted in Fig. 3.

3.2 Experimental phantom results

3.2.1 Layered phantom results

The predicted results of the CNN and XU are the average across the collected frames. The uncertainty is measured by dividing the standard deviation of the results by their mean value. The results of the patchless CNN and XU algorithm are given in Fig. 4. The uncertainty map has a higher value in the region where the CNN incorrectly estimates a lower value in the phantom (the region is marked with red arrows).

In order to quantitatively investigate the presented results, the ratio of the mean of obtained values in regions R2 and R1 (highlighted in the B-mode image) are computed and compared with the ratio of their known backscattering coefficient. The ratio is $\frac{6.37}{3.52}=1.81$ for the backscattering coefficient, it is $\frac{16.34\pm 3.72}{10.99\pm 2.73}=1.61\pm 0.63$ for XU, and $\frac{5.19\pm 0.31}{1.81\pm 0.03}=1.81\pm 0.11$ for the patchless CNN. It can be seen that the ratio for the patchless CNN is closer to the true backscattering coefficient ratio compared to the XU method which demonstrates that this method has a lower bias compared to XU.

3.2.2 CIRS phantom results

The results are shown in Fig. 5. The XU algorithm fails to obtain reliable values on the borders (marked in b) of the inclusions due to presence of two different distributions inside the patches of those regions, while the CNN results do not suffer from this issue. The CNN provides higher contrast and less bias with the samples with lower correlation, but due to skipping more samples, the inclusion borders are not as clear as with the higher correlation one. This demonstrates a trade off between skipping the samples to reduce the correlation (lower bias) and spatial resolution. This issue was less identifiable in the layered phantom since there were adequate samples in each layer; therefore, we were able to skip enough samples to reduce the bias.

3.2.3 Gammex phantom results

The results are shown in Fig. 6. The parametric image obtained by XU algorithm has artifacts around the inclusions similar to Fig. 5. The parametric image with higher correlation (c) has higher spatial resolution but it has more biased values (compare the inclusion 1 in (c) and (d)). The inclusion 1 is clearly detected by the CNN for both inputs (c and d), but in (c) the inclusion has a more biased values than (d) due to presence of higher correlation. The inclusion 2 is also hardly detected by the CNN since only part of it is available.

4 Discussion

In this paper, a deep learning frame work is proposed to quantify the scatterer number density. One important thing to note is that the method of calculating the scatterer number density in simulation data generation dictates the output of the network. In other word, different definition of resolution cell in simulation data generation can change the ground truth value which also alters the output of the network.

Another crucial point is that the correlation between samples can affect the results. Most of recent work focus on estimating the scatterer number density using training data obtained from i.i.d samples. Although these samples contain more information than the correlated ones, achieving this condition in real data is difficult and requires skipping many samples which is not efficient in parametric image reconstruction. In this work, the approach is changed by generating simulation data using the grid-based method to have correlation between the samples. The disadvantage of this approach is that the exact value of $\alpha$ is not known, and generating coherent component is complex.

The impact of different correlations for the experimental phantoms is also investigated. Generally, lower correlation can help having less bias in estimation but achieving this might cause losing spatial resolution due to requiring skipping many samples.

4.1 Conclusion

In this paper, a patchless deep learning solution was proposed to obtain scatterer number density parametric images. The method was validated using simulation data and three experimental phantoms. Frame-wise uncertainty map was also obtained from the parametric images of the multiple frames. In order to enable the network to predict the experimental phantom data having correlated samples, a grid-based simulation in which correlated samples were present was employed for training data generation. The proposed method was trained on simulation data and was able to reconstruct scatterer number density parametric image of experimental phantoms imaged by different scanners.

4.2 Acknowledgments

We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), and the report of the layered phantom data in the previous publication [9].