---
abstract: |
    Removing speckle noise from medical ultrasound images while preserving image features without introducing artifact and distortion is a major challenge in ultrasound image restoration. In this paper, we propose a multiframe-based adaptive despeckling (MADS) algorithm to reconstruct a high-resolution B-mode image from raw radio-frequency (RF) data that is based on a multiple input single output (MISO) model. It utilizes the speckle patterns estimated using a novel multiframe-based adaptive approach for ultrasonic speckle noise estimation (MSNE) based on a single input multiple output (SIMO) modeling of consecutive deconvolved ultrasound image frames. The elegance of the proposed despeckling algorithm is that it addresses the despeckling problem by completely following the signal generation model unlike conventional ad-hoc smoothening or filtering based approaches, and therefore, it is likely to maximally preserve the image features. As deconvolution is a necessary pre-processing step to despeckling, we describe here a $2$-D extension of the SIMO model-based $1$-D deconvolution method. Finally, a complete framework for the generation of high-resolution ultrasound B-mode image has been also established in this paper. The results show $8.55-15.91$ dB, $8.24-14.94$ dB improvement in terms of SNR and PSNR, respectively, for simulation data and $2.22-3.17$ improvement in terms of NIQE for *in-vivo* data compared to the traditional despeckling algorithms. Visual comparison shows superior texture, resolution, details of B-mode images offered by our method compared to those by a commercial scanner, and hence, it may significantly improve the diagnostic quality of ultrasound images.

    Clinical ultrasound, analytic modeling, SIMO model, $2$-D deconvolution, MISO model, despeckling, high-resolution B-mode.
author:
- 'Jayanta Dey, Sharmin R. Ara, and Md. Kamrul Hasan'
bibliography:
- 'jd.bib'
title: |
    Multiframe-based Adaptive Despeckling Algorithm for Ultrasound B-mode Imaging with Superior Edge and Texture

    [^1]
---

INTRODUCTION
============

(US) imaging system being non-invasive, non-ionizing, portable, and cost effective has become the most prevalent diagnostic tool among all the currently available imaging modalities, e.g., X-ray, magnetic resonance imaging, and computed tomography. However, imperfection of US imaging system design and the underlying physical phenomena related with US image acquisition give rise to low resolution and speckle noise that tend to reduce the image contrast, obscure and blur image details such as inclusion and small structure boundary, tissue texture and thereby, decrease the quality and reliability of medical ultrasound [@michailovich2006despeckling]. Speckle noise is a granular pattern inherent in any coherent imaging modalities [@choi2018speckle], [@de2011real] similar to ultrasound imaging. It results from the constructive and destructive interferences of the reflected echos with different phases and amplitudes from the target at the receiving transducer. Removal of speckle noise from US images is difficult due to its multiplicative nature and the challenge of maintaining the precise texture of the image [@joel2018extensive]. As traditional despeckling filters distort the original image texture and introduce artifacts like blurring edges, changing the shape of structures present in the image by smoothening the noise corrupted image, the original noise affected images are sometimes more preferred than the noise-removed ones in the analysis where the image details have high importance [@joel2018extensive]. Therefore, an effective signal processing approach for removing speckle noise from the image while preserving the original tissue texture and small details of the image is vital to increase the diagnostic potential of medical ultrasound.

A number of algorithms have been reported in the literature for despeckling US images which attempts to remove speckle noise either in their mulplicative form or converting the noise in the additive form. Among the first group of algorithms, the spatial averaging based approaches exploit the repetitive nature of the US image, and among them, linear filtering methods, such as Gaussian filter and mean filter are effective in reducing speckle noise [@choi2018speckle]. However, they tend to oversmooth the texture and blur edges present in the image. To overcome this problem, nonlinear approaches based on local [@frost1982model], [@kuan1985adaptive] and non-local statistics [@coupe2009nonlocal], [@buades2005review] of the image have been proposed. These algorithms are mainly weighted filters, in which the weights depend on the similarity between the intensity values of the patches surrounding the pixels [@shahdoosti2018maximum]. The main difference between local and non-local means methods is that the non-local means method employs the most similar pixels in the image to denoise the current pixel regardless of their Euclidean distance. Although these approaches tend to preserve textures and edges, their performance is dependent on tuning parameters, such as filter and patch size. Again, the nonlinear approaches based on the diffusion equation [@yu2002speckle] not only preserves edges but also enhances edges by inhibiting diffusion across edges and allowing diffusion on either side of the edge. However, selection of the parameter-values is a major issue in this method, as a value of parameter that is smaller than the optimum one leads to unsatisfactory noise suppression whereas a higher parameter-value results in poor structure preservation [@jain2019non].

The second group of algorithms converts the multiplicative speckle noise into an additive one by using logarithmic transformation. Now, as shown in [@michailovich2006despeckling], the additive noise can be handled using any traditional denoising scheme, and the performance of this scheme determines the overall efficacy of despeckling. Finally, the denoised image is exponentially transformed back to give the despeckled image. The overall process is termed as homomorphic filtering approach [@michailovich2006despeckling], and among the denoising schemes used in this process, the transfrom domain or multi-resolutional based appraoches [@gupta2005locally; @gupta2005despeckling; @yue2006nonlinear; @koundal2016speckle] are of higher efficacy. In [@luisier2007new], an advanced ultrasound despeckling algorithm is proposed based on the intra-scale correlation between the wavelet coefficients. Among the multi-resolutional approaches, as shown in [@joel2018extensive], despeckling based on non-subsampled contourlet gives superior performance. However, in wavelet-based despeckling, a threshold is a critical parameter that is to be determined based on the *a priori* knowledge of the distribution of the speckle pattern. In addition, the threshold-based filtering of wavelet coefficients implies texture smoothening, and it gives rise to artifact such as Gibbs phenomenon near the edges [@bhuiyan2007spatially]. All of the approaches discussed so far relies on ad-hoc filtering or smoothening technique without addressing mathematically the speckle noise generation model. Therefore, it cannot be guaranteed that these algorithms only operate on the speckle noise without significantly distorting the true image. In this paper, we propose a multiframe-based adaptive despeckling (MADS) algorithm that treats the speckle noise in its multiplicative form and utilizes the speckle patterns estimated using the multiframe-based adaptive ultrasonic speckle noise estimation (MSNE) algorithm proposed in this paper. The MSNE algorithm is based on formulating the true image as single input and the envelope of deconvolved consecutive US image frames with multiplicative speckle noise pattern in each frame as multiple outputs. The despeckling algorithm, on the contrary, treats the envelope of deconvolved consecutive US image frames as multiple inputs and the true image as single output. According to the mathematical model representing the US imaging system, deconvolution is necessary prior to despeckling for resolution enhancement of the raw RF data. Hence, a $2$-D deconvolution approach as an extension of our previously published $1$-D deconvolution algorithm, i.e., *b*MCFLMS [@dey2018ultrasonic] has been also described in this paper. To prevent misconvergence of the MSNE algorithm in the presence of additive noise and estimation error resulting from the deconvolution step, a zero-lag correlation contraint derived from the deconvolved image and the estimated speckle pattern is attached with the original cost function. As the overall despeckling approach has been derived by completely following the signal generation models, it is likely to maximally preserve the diagnostically important details and tissue texture present in the image. Finally, a complete framework including deconvolution, despeckling, and post-processing for ultrasound B-mode image generation with superior quality in terms of resolution, edges, small details and texture is also established in this paper for greater interest of the researchers.

The paper is organized as follows. Section II describes the US imaging model and formulates the problems to be addressed. The $2$-D deconvolution and speckle noise estimation based on SIMO models along with the derivation of a MISO model for despeckling the decovolved US images using the estimated speckle patterns are presented in Section III. Post-processing on the despeckled image is explained in Section IV. The performance of the proposed method is tested using simulation and *in-vivo* data in Section V, and the pros and cons of the proposed framework is discussed in Section VI. Finally, summarizing the contributions with highlights for future research the paper concludes in Section VII.

Problem Formulation
===================

Low resolution and speckle noise are the major issues related with US imaging. However, an overall US image enhancement can be achieved by addressing these two issues in a sequential two-step process as described in [@michailovich2006despeckling]. First, the correlation between the image samples is to be minimized to increase the image resolution, and second, speckle noise has to be removed from the decorrelated image to improve the image contrast and better visualize the tissue texture. To achieve the aforementioned two objectives, a suitable model representing the US imaging system is necessary. With this in view, considering linear wave propagation through the tissue, and the scattering of the ultrasound pulse in the tissue as weak, we can use the first order Born approximation and consider the tissue scattering system as a linear system [@michailovich2006despeckling]. Therefore, the blurring, i.e., low resolution of an RF-image can be modeled as the result of convolution between the point-spread function (PSF) $s(m,n)$ of the imaging system with the tissue reflectivity function (TRF) $h(m,n)$ [@michailovich2006despeckling], [@xue2018parametric], [@al2018faster]. Mathematically, this can be written as $$\label{convolutionModel}
x(m,n) = s(m,n)*h(m,n) + v(m,n)$$ where $x(m,n)$ is the backscattered ultrasound image data from the $n$-th A-line at discrete time $m$ and $v(m,n)$ is the additive noise associated with measurement error and other physical phenomena not accounted by the convolution model. In US imaging, $h(m,n)$ is corrupted by speckle noise, and as described in [@michailovich2006despeckling], [@zong1998speckle], [@achim2001novel], $h_e(m,n)$, i.e., the envelope of $h(m,n)$ can be modeled as multiplicative with the true image as $$\label{speckleModel}
h_e(m,n) = r(m,n) u(m,n) + \zeta(m,n)$$ where $r(m,n)$, $u(m,n)$, and $\zeta(m,n)$ are true image, speckle noise, and additive noise resulting from the deconvolution process and the portion not accounted by the multiplicative process, respectively. Therefore, in ultrasound imaging, we are given the back-scattered data $x(m,n)$, and our objective is to design algorithms to estimate the TRF $h(m,n)$ from $x(m,n)$, and thereafter, obtain the despeckled image $r(m,n)$ from the envelope of $h(m,n)$.

Method
======

In this paper, our main concern is to derive a novel despeckling algorithm that is likely to preserve maximum features present in the image. However, according to and , deconvolution is a necessary pre-processing step to despeckling. Therefore, for the completeness of a high-resolution B-mode image generation, we also consider it important to include a $2$-D extension of our previously proposed $1$-D deconvolution (*b*MCFLMS) [@dey2018ultrasonic] algorithm and some post-processing techniques like gamma correction and grey level transformation so that the whole paradigm is described in a single paper.

Deconvolution of RF Echo Data
-----------------------------

In this section, we attempt to estimate the TRF, $h(m,n)$, in with increased resolution by removing the effect of the PSF, $s(m,n)$, from the raw RF data, $x(m,n)$.

However, to estimate $h(m,n)$, similar to the approach as described in [@abeyratne1995higher], we consider the $2$-D distortion kernal, i.e., PSF $s(m,n)$ in decomposible into two $1$-D distortion kernals (PSFs): one along the axial direction and the other along the lateral direction. Following this assumption, can be modified as $$x(m,n) = s_a(m) *_a s_l(n) *_l h(m,n) + v(m,n)$$ where $s_l(n)$ and $s_a(m)$ are the lateral and axial PSFs, respectively, and ‘$*_l$’ and ‘$*_a$’ represent convolution along the lateral and the axial directions, respectively. A novel technique for removing the effect of axial PSF $s_a(m)$ from the measured RF image $x(m,n)$ was reported in our previous work [@dey2018ultrasonic] using a single input multiple output (SIMO) model as shown in Fig. \[deconvolution\], where $s_a(m)$ convolves with the $i$-th A-line of the axial TRF denoted as $h_i^a(m)$ and with additive noise $v_i(m)$ gives the $i$-th A-line RF data $x_i(m)$: $$\label{axialConv}
x_i(m) = s_a(m) *_a h_i^a(m) + v_i(m)$$ where $$h_i^a(m) = s_l(n) *_l h(m,n)$$ In matrix form, can be written as $$\mathbf{X} = \mathbf{S}_a \mathbf{h}^a + \mathbf{v}$$ where $\mathbf{S}_a$ is the convolution matrix formed using the axial PSF $s_a(m)$ and $$\nonumber
\mathbf{X} = 
\begin{bmatrix}
\mathbf{X}_1 & \mathbf{X}_2 & \cdots & \mathbf{X}_{N'}
\end{bmatrix},$$ $$\nonumber
\mathbf{h}^a = 
\begin{bmatrix}
\mathbf{h}_1^a & \mathbf{h}_2^a & \cdots & \mathbf{h}_{N'}^a
\end{bmatrix}.$$ Here, $N'$ is the total number of A-lines, and $\mathbf{X}_i$ and $\mathbf{h}_i^a$ are the $i$-th A-line with $M'$ samples taken from $x_i(m)$ and $h_i^a(m)$, respectively where $m = 1, 2, \cdots, M'$. To account for the non-stationarity of the axial PSF, the RF data were divided into $B$ blocks with equal length $L_b$, and a block-based cost function $J^b$ for the $b$-th block was formulated (for details see [@dey2018ultrasonic]) to estimate the axial TRF block-by-block in the frequency-domain as $$J^b = \sum_{i=1}^{n-1} \sum_{j=i+1}^n \tilde{\underline{\mathbf{e}}}_{ij}^{bH} \tilde{\underline{\mathbf{e}}}_{ij}^{b}$$ where, ‘$H$’ denotes the Hermitian operation, any variable with ‘$\underline{~~}$’ represents the variable in the frequency-domain, and $\tilde{\underline{\mathbf{e}}}_{ij}^{b}$ is the Fourier transform of $\tilde{\mathbf{e}}_{ij}^{b}$ defined as $$\nonumber
\tilde{\mathbf{e}}_{ij}^b = \sum_{p=1}^{b-1} \mathbf{A}_1 \mathbf{e}_{ij}^{p(b-p)} + \sum_{p=1}^b \mathbf{A}_2 \mathbf{e}_{ij}^{p(b-p+1)}$$ Here, $$\label{eqn:last_err}\nonumber
\mathbf{e}_{ij}^{p(b-p)} = \mathbf{C}_{\tilde{x}_i^p}\hat{\mathbf{h}}_j^{a(b-p)} - \mathbf{C}_{\tilde{x}_j^p}\hat{\mathbf{h}}_i^{a(b-p)} ,
p = 1,2,\cdots,b-1,$$ $$\label{eqn:first_err}\nonumber
\mathbf{e}_{ij}^{p(b-p+1)} = \mathbf{C}_{\tilde{x}_i^p}\hat{\mathbf{h}}_j^{a(b-p+1)} - \mathbf{C}_{\tilde{x}_j^p}\hat{\mathbf{h}}_i^{a(b-p+1)},
p = 1,2,\cdots,b,$$ and $\mathbf{A}_1$, $\mathbf{A}_2$ are the truncation matrices truncating the last $(L_b-1)$ and the first $L_b$ samples of the error function, respectively. $\mathbf{C}_{\tilde{x}_i^p}$ is the convolution matrix formed using the RF data, $x(m,n)$ along the $i$-th A-line and the $p$-th block. Now, the $b$-th block axial TRF $\underline{\mathbf{h}}^{ab}$ was estimated as $$\label{eqn:arg}
\underline{\hat{\mathbf{h}}}^{ab}  = arg_{\underline{\mathbf{h}}^{ab}} \min J^b, \mbox{subject~to~}
||\underline{\hat{\mathbf{h}}}^a|| = 1$$ where ‘$||\cdot||$’ denotes the $l_2$-norm and $$\label{eqn:total_trf}
\underline{\hat{\mathbf{h}}}^a =
\begin{bmatrix}
\underline{\hat{\mathbf{h}}}^{a1T} & \underline{\hat{\mathbf{h}}}^{a2T} & \cdots  & \underline{\hat{\mathbf{h}}}^{aBT}
\end{bmatrix}^T$$ with ‘$T$’ denoting matrix transpose operation and $$\underline{\hat{\mathbf{h}}}^{ab} =
\begin{bmatrix}
\underline{\hat{\mathbf{h}}}_{1}^{ab} & \underline{\hat{\mathbf{h}}}_{2}^{ab}  & \cdots  & \underline{\hat{\mathbf{h}}}_{N'}^{ab}
\end{bmatrix}$$ Here, $\underline{\hat{\mathbf{h}}}_{i}^{ab}$ denotes the estimated $i$-th A-line $b$-th block axial TRF. In sample-domain, the estimated axial TRF along the $i$-th A-line can be written as $\hat{h}_i^a(m)$. Now, the estimated axial TRF $\hat{h}_i(m)$ along the $i$-th A-line at discrete time $m$ can be modeled as $$\label{lateralDeconv}
\hat{h}_i^a(m) = s_l(n) *_l h(m,n) + v'(m,n)$$ where $v'(m,n)$ is the noise resulting from the estimation error of the axial TRF. Therefore, from and , it is apparent that an attempt, similar to the approach adopted in the axial direction, can be made in the lateral direction to undo the effect of the lateral distortion kernal (PSF) from the estimated axial TRF $\hat{h}_i^a(m)$ to estimate the $2$-D deconvolved TRF $h(m,n)$. In this approach, the PSF is considered to be laterally stationary as described in [@michailovich2005novel], and therefore, no blocking is required in the lateral direction. The method is summarized in Table \[two-D-mcflms\].

------------------------------------------------------------------------

1.  1.  Select the *b*MCFLMS method reported in [@dey2018ultrasonic]

2.  1.  Set $\mathbf{X} =$ raw radio-frequency (RF) data, where $\mathbf{X}$ is the data to be deconvolved

    2.  Set the data length, $L = $ axial length, and the number of channels, $M = $ lateral length of the raw RF data. Execute the *b*MCFLMS algorithm with block number, $B = 2$ along the axial direction.

3.  1.  Set $\mathbf{X} =$ axially deconvolved data in step 2

    2.  Set the data length, $L = $ lateral length, and the number of channels, $M = $ axial length of the raw radio-frequency (RF) data. Execute the *b*MCFLMS algorithm with block number, $B = 1$ along the lateral direction.

------------------------------------------------------------------------

The estimated TRF after lateral deconvolution is given by $$\underline{\hat{\mathbf{h}}} =
\begin{bmatrix}
\underline{\hat{\mathbf{h}}}_1 & \underline{\hat{\mathbf{h}}}_2 & \cdots  & \underline{\hat{\mathbf{h}}}_{M'}
\end{bmatrix}^T$$ where, $\underline{\hat{\mathbf{h}}}_i$ is the estimated lateral TRF along the $i$-th sample line, i.e., samples along the $i$-th row of $\underline{\hat{\mathbf{h}}}$. In sample-domain, the $m$-th row and $n$-th column sample of the estimated TRF, $\underline{\hat{\mathbf{h}}}$ can be denoted as $\hat{h}(m,n)$.

Proposed Despeckling Algorithm {#sssec:ref1}
------------------------------

The envelope of the estimated TRF $\hat{h}(m,n)$ in the previous subsection is corrupted with speckle noise as given by . It is apparent from that knowledge of the speckle pattern $u(m,n)$ in an image frame can help despeckling that frame. However, in the absence of additive noise, direct division of $h_e(m,n)$ by $u(m,n)$ may amplify noise and/or give rise to division by zero problem. Hence, in this subsection, we attempt to formulate an energy constrained iterative approach to find an speckle noise cancelling multiplying factor to despeckle the frame. In what follows, we attempt to formulate a novel SIMO model for the deconvolved image frames and thereby, estimate the speckle noise in the respective deconvolved frames using an adaptive filtering technique. Then a novel MISO model will be proposed to despeckle the decconvolved image frames using the estimated speckle pattern of the respective frames.

In an ultrasound imaging system, images are generally acquired at a frame rate ranging from $10-60$ frames per second (FPS) with speckle patterns generated randomly in each image frame from the interference of the US pulse at the receiving transducer. In our work, we attempt to use $p$ consecutive deconvolved frames to formulate a single input multiple output (SIMO) model as shown in Fig. \[speckle\]. Here, the true ultrasound image $r(m,n)$ that is considered stationary throughout the $p$ frames, multiplies with the speckle noise $u_i(m,n)$ of the $i$-th frame, and with an additive noise $w_i(m,n)$ gives the envelope, $\hat{h}_{ei}(m,n)$, of the estimated deconvolved image $\hat{h}_i(m,n)$ of the $i$-th frame: $$\label{simoBasic}
\hat{h}_{ei}(m,n) = r(m,n) u_i(m,n) + w_i(m,n)$$ In matrix form, can written as

$$\hat{\mathbf{H}}_{ei} = \mathbf{R}\cdot*\mathbf{U}_i + \mathbf{W}_i$$

where ‘$\cdot*$’ denotes elementwise multiplication. Here, $\hat{h}_{ei}(m,n)$, $r(m,n)$, $u_i(m,n)$ and $w_i(m,n)$ represents the $m$-th row and $n$-th column elements of $\hat{\mathbf{H}}_{ei}$, $\mathbf{R}$, $\mathbf{U}_i$ and $\mathbf{W}_i$, respectively. Therefore, here the challenge is to estimate the speckle noise $u_i(m,n)$ from each frame in the presence of additive noise $w_i(m,n)$ and then remove the speckle noise from the deconvolved envelope image $\hat{h}_{ei}(m,n)$. The assumptions behind the SIMO model formulation and the identifiability condition [@xu1995least] for the speckle noise pattern $u_i(m,n)$ of the $i$-th frame are

1.  The true image $r(m,n)$ is stationary throughout the $p$ consecutive frames.

2.  The speckle patterns of each frame do not share common zeroes with the rest $p-1$ consecutive frames.

These assumptions are realistic because for an ultrasound video recording with $30$ FPS, consecutive $5-10$ frames take around $0.17-0.33~$second during which the hand motion can be ignored. Then we can consider the true ultrasound image $r(m,n)$ as stationary throughout these frames. Again, formation of speckle noise in consecutive frames is a completely random phenomena, and hence, they are unlikely to contain common zeros. However, the probablity of sharing common zeros between the speckle patterns of the consecutive frames can be further reduced by increasing the number of fames into consideration; doing this will improve the identification accuracy of the patterns provided that the stationarity assumption remains valid as shown in the result section later. In what follows, we derive a multiframe-based adaptive speckle noise estimation algorithm using the proposed SIMO model.

### Speckle Noise Estimation

In the absence of additive noise, the following error function $e_{ij}(m,n)$ can be used to estimate the speckle noise: $$\label{error}
e_{ij}(m,n) = \hat{h}_{ei}(m,n) \hat{u}_j(m,n) - \hat{h}_{ej}(m,n) \hat{u}_i(m,n)$$ where $\hat{u}_i(m,n)$ is the estimated speckle noise of the $i$-th frame. Notice that for additive noiseless case if we can estimate the speckle pattern accurately, the error function defined in becomes zero. Using this fact, we can build the following cost function to iteratively estimate the speckle noise:

$$\label{costFunction}
J = \sum_{i=1}^{p-1} \sum_{j = i+1}^p ||\mathbf{E}_{ij}\cdot* \mathbf{E}_{ij}||_F^2$$

where $$\mathbf{E}_{ij} = \hat{\mathbf{H}}_{ei}\cdot* \hat{\mathbf{U}}_j - \hat{\mathbf{H}}_{ej}\cdot*\hat{\mathbf{U}}_i$$ and ‘$||\cdot||_F$’ indicates the Frobenius norm. An estimate of the speckle noise $\hat{\mathbf{U}}$ can be obtained by minimizing the cost function $J$ as

$$\label{estimateU}
\hat{\mathbf{U}}  = arg_{\mathbf{U}} \min J, \mbox{subject~to~}||\hat{\mathbf{U}}||_F = 1$$

where $$\label{s}
\hat{\mathbf{U}} =
\begin{bmatrix}
\hat{\mathbf{U}}_1 & \hat{\mathbf{U}}_2  & \cdots  & \hat{\mathbf{U}}_p
\end{bmatrix}$$ Taking the gradient of $J$ in , we get $$\label{gradient}
\nabla_k J = \frac{\partial J}{\partial \hat{\mathbf{U}}_k} = 2\sum_{i=1}^p \hat{\mathbf{H}}_{ei}.* \mathbf{E}_{ik}$$ The update equation for the MSNE algorithm at the $q$-th iteration is

$$\label{update}
\hat{\mathbf{U}}(q+1) = \hat{\mathbf{U}}(q) - \mu(q) \nabla J(q)\big|_{\mathbf{U}=\hat{\mathbf{U}}(q)}$$

where, $$\begin{aligned}
\nabla J(q) &= \frac{\partial J(q)}{\partial \hat{\mathbf{U}}} \nonumber\\
&=
\begin{bmatrix}
\nabla_1 J(q) & \nabla_2 J(q)  & \cdots  & \nabla_p J(q)
\end{bmatrix}\end{aligned}$$ and $\mu(q)$ is the variable step-size (VSS) which is such that the misalignment of $\hat{\mathbf{U}}(q+1)$ with the true noise pattern $\mathbf{U}$ is minimum at every iteration, given the current estimate $\hat{\mathbf{U}}(q)$: $$\begin{aligned}
\label{muUp}
J_\mu (q) &= (|| \mathbf{U} - \alpha \hat{\mathbf{U}}(q+1) ||_F)^2| \hat{\mathbf{U}}(q) \nonumber\\
          &= \big|\big| \mathbf{U}\cdot*\mathbf{U} - 2\alpha \mathbf{U}\cdot*\hat{\mathbf{U}}(q) + 2\alpha\mu (q)\mathbf{U}\cdot*\nabla J(q) \nonumber \\
          & + \alpha^2\hat{\mathbf{U}}(q)\cdot*\hat{\mathbf{U}}(q) - 2\alpha^2\mu (q)\hat{\mathbf{U}}\cdot*\nabla J(q) \nonumber\\
          & + \alpha^2\mu^2 (q) \nabla J(q)\cdot*\nabla J(q) \big|\big|_S | \hat{\mathbf{U}}(q)\end{aligned}$$ where $\alpha$ is a scaling constant inherent in any blind channel identification approach based on the cross-relation, and we define an operator ‘$||\cdot||_S$’ which evaluates the sum of the matrix elements. Minimizing , i.e., setting the gradient of $J_\mu(q)$ with respect to $\mu(q)$ to zero, we get $\mu(q)$ as $$\label{vss}
\mu(q) = \frac{\big|\big|\hat{\mathbf{U}}(q)\cdot*\nabla J(q)\big|\big|_S}{\big|\big|\nabla J(q)\cdot*\nabla J(q)\big|\big|_S}$$ Equation can be considered as a variant of VSS derived in [@gaubitch2006generalized].

Thus far additive noise has been ignored in the derivation of the proposed MSNE algorithm. However, it has a similar effect on the convergence of the MSNE algorithm as described in [@dey2018ultrasonic]. To solve the problem, we need to impose a constraint on so as to prevent the deviation of the estimated speckle pattern from the true speckle pattern. To this end, consider the following model of speckle corrupted image for the $k$-th frame as described in [@coupe2009nonlocal]: $$\label{simEqn}
h_{ek}(m,n) = r(m,n) + r^\eta(m,n)~\mathcal{V}_k(m,n)$$ where $\mathcal{V}_k(m,n) \sim \mathcal{N}(0,~\sigma^2)$. Ignoring the additive noise in and comparing with , our proposed algorithm is basically estimating the true speckle pattern as $$\label{trueSpeckle}
u_k(m,n) = 1 + r^{(\eta-1)}(m,n)~\mathcal{V}_k(m,n)$$ Then the speckle pattern $\hat{u}_k(m,n)$, estimated using the proposed method, can be expresses as $$\label{estimatedSpeckle}
\hat{u}_k(m,n) = 1 + r^{(\eta-1)}(m,n)~\hat{\mathcal{V}}_k(m,n)$$ Since $r(m,n)$ is not a variable here, to make the estimated speckle pattern $\hat{u}_k(m,n)$ in close to the true speckle pattern $u_k(m,n)$ in , we need to maximize the zero-lag correlation between $\hat{\mathcal{V}}_k(m,n)$ and $\mathcal{V}_k(m,n)$. However, to attach this as a constraint on , estimates of the true image $r(m,n)$ and the parameter $\eta$ are necessary. Alternatively, consider the zero-lag correlation between the deconvolved image $\hat{h}_{ek}(m,n)$ and the estimated speckle pattern $\hat{u}_k(m,n)$: $$\begin{aligned}
\label{constraint}
J_{corr} &= \big| \big| \hat{\mathbf{H}}_{ek}\cdot*\hat{\mathbf{U}}_k \big| \big|_S \nonumber \\
         &= \big| \big| (\mathbf{H}_{ek} + \mathbf{W}_k)\cdot*\hat{\mathbf{U}}_k \big| \big|_S \nonumber \\
         &= \big| \big| \mathbf{H}_{ek}\cdot*\hat{\mathbf{U}}_k \big| \big|_S\end{aligned}$$ where the zero-lag correlation between the additive noise $\mathbf{W}_k$ and speckle noise $\mathbf{U}_k$ is considered zero. Using and in , we get $$\begin{aligned}
\label{show}
J_{corr} &= \big| \big| \mathbf{R} + \mathbf{R}^\eta\cdot*\mathcal{V}_k + \mathbf{R}^\eta\cdot*\hat{\mathcal{V}}_k + \mathcal{V}_k\cdot* \hat{\mathcal{V}}_k \big| \big|_S \nonumber\\
         &= \big| \big| \mathbf{H}_{ek} + \mathbf{R}^\eta\cdot*\hat{\mathcal{V}}_k + \mathcal{V}_k\cdot* \hat{\mathcal{V}}_k \big| \big|_S \nonumber\\
         &= c + \big| \big| \mathcal{V}_k\cdot* \hat{\mathcal{V}}_k \big| \big|_S\end{aligned}$$ where at the point of misconvergence when the estimated speckle pattern is close to the true speckle pattern, the zero-lag correlation between $r^\eta(m,n)$ and $\hat{\mathcal{V}}_k(m,n)$ can be considered zero, and $c$ is a constant defined as $c = \big| \big| \mathbf{H}_{ek} \big| \big|_S$. From , it is apparent that maximizing $J_{corr}$ or equivalently minimizing $-J_{corr}$ is analogous to maximizing the zero-lag correlation between $\hat{\mathcal{V}}_k(m,n)$ and $\mathcal{V}_k(m,n)$. To prevent the misconvergence of the proposed algorithm in the noisy case, we propose to use the zero-lag correlation constraint in with the MSNE cost function in . Then modifying for noisy case, we obtain $$\label{constraint2}
J_c(q) = J(q) - \beta_1 J_{corr}(q)$$ where $\beta_1$ is the Lagrange multiplier or also known as the coupling factor. Taking gradient of with respect to $\hat{\mathbf{U}}_k$, we get $$\label{gradientConstraint}
\nabla_k J_c(q) = \nabla_k J(q) - \beta_1 \mathbf{H}_{ek}$$ Replacing $\nabla_k J(q)$ in by $\nabla_k J_c(q)$, we can estimate the speckle noise pattern in each of the $p$ frames.

### Estimation of the True Ultrasound Image

Hitherto, a novel algorithm for estimating the speckle pattern has been explained with a view to estimating the true ultrasound image $r(m,n)$ in using $\hat{\mathbf{U}}$ from . Now, we describe a novel MISO model as shown in Fig. \[m-MINT\], where the speckle noise cancellation (SNC) factors for each of the $i$-th frame $g_i(m,n)$ multiplies with $h_{ei}(m,n)$, to obtain an estimate of the true US image $\hat{r}_i(m,n)$ for that frame. In the absense of additive noise in and estimation error in the estimated speckle pattern $\hat{\mathbf{U}}_i$, the estimated ultrasound image of the $i$-th frame, $\hat{\mathbf{R}}_i$, can be obtained in matrix form as

$$\begin{aligned}
\hat{\mathbf{R}}_i &= \mathbf{H}_{ei}\cdot*\mathbf{G}_i \label{equalize}\\
                    &= \mathbf{R}_i\cdot*\mathbf{U}_i\cdot*\mathbf{G}_i \label{simpleEqualize}\end{aligned}$$

where $\mathbf{G}_i$ is the elementwise multiplying SNC factor to equalize the speckle pattern for the $i$-th frame. Here for an accurate estimation of $\mathbf{R}_i$, the elements of the matrix $\mathbf{U}_i.*\mathbf{G}_i$ in should be equal to $1$. Therefore, we can formulate the following cost function to estimate the SNC factors: $$\label{MINTcost}
J_{eq} = ||\mathbf{U}\cdot*\mathbf{G} - \mathbf{D}||_F^2$$ where $$\label{g}
\mathbf{G} =
\begin{bmatrix}
\mathbf{G}_1 & \mathbf{G}_2  & \cdots  & \mathbf{G}_p
\end{bmatrix}$$ and $\mathbf{D}$ is a matrix with all entries equal to $1$. Taking gradient of with respect to $\mathbf{G}$, we get $$\label{gradMINT}
\nabla J_{eq} = \frac{\partial J_{eq}}{\partial \mathbf{G}} = 2(\mathbf{U}.*\mathbf{G} - \mathbf{D}).*\mathbf{U}$$ The update equation for estimating $\mathbf{G}$ at the $q'$-th iteration is given by $$\label{updateMINT}
\mathbf{G}(q'+1) = \mathbf{G}(q') - \mu(q') \nabla J_{eq}(q')$$ where $\mu(q')$ is the VSS for the $q'$-th iteration, and can be obtained following as $$\label{mu}
\mu(q') = \frac{\big|\big|\mathbf{G}(q')\cdot*\nabla J_{eq}(q')\big|\big|_{S}}{\big|\big|~\nabla J_{eq}(q')\cdot*\nabla J_{eq}(q')\big|\big|_S}$$ Up to now, we have ignored the effect of estimation error in $\hat{\mathbf{U}}$. If we assume that $\mathbf{E}_s$ be the estimation error in $\hat{\mathbf{U}}$, then we can write $$\label{errorU}
\hat{\mathbf{U}} = \mathbf{U} + \mathbf{E}_s$$ Now, replacing the true speckle pattern $\mathbf{U}$ with the estimated speckle pattern $\hat{\mathbf{U}}$, the cost function in becomes $$\begin{aligned}
J_{eq} &= ||\hat{\mathbf{U}}.*\mathbf{G} - \mathbf{D}||_F^2 \label{hatEqualize}\\
       &= ||(\mathbf{U}+\mathbf{E}_s).*\mathbf{G} - \mathbf{D}||_F^2 \nonumber\\
       &= ||\mathbf{U}.*\mathbf{G} - \mathbf{D}||_F^2 + \mathcal{E}\end{aligned}$$ where $\mathcal{E}$ represents the terms including the estimation error $\mathbf{E}_s$. Therefore, the gradient in will have two components– one from the desired part of the cost function and the other from the estimation error part, i.e., $$\label{gradMINT2}
\nabla J_{eq}(q') = \nabla J_{eq}^{desired}(q') + \nabla J_{eq}^{error}(q')$$ From , we observe that $\hat{\mathbf{U}}$ and $\mathbf{G}$ have inverse relation so as to make their elementwise product equal to $1$. Therefore, in , the term associated with the estimation error in $\hat{\mathbf{U}}$ gives rise to an SNC factor $\mathbf{G}$ in which the small estimation error is magnified. To make $\mathbf{G}$ less sensitive to such phenomenon, we impose an energy regularization constraint on $\mathbf{G}$ in : $$\label{modifiedCost}
J'_{eq} = ||\hat{\mathbf{U}}\cdot*\mathbf{G} - \mathbf{D}||_F^2 + \beta_2||\mathbf{G}||_F^2$$ where $\beta_2$ is the Lagrange multiplier. Now, the gradient in becomes $$\label{gradMINT3}
\nabla J'_{eq} = 2(\hat{\mathbf{U}}\cdot*\mathbf{G} - \mathbf{D})\cdot*\hat{\mathbf{U}} + 2 \beta_2\mathbf{G}$$ In addition to the energy regularization constraint, we propose the following gradient averaging technique to average out or at least reduce the detrimental effect of $\nabla J_{eq}^{error}(q')$: $$\nabla J''_{eq}(q') = \alpha \nabla J'_{eq}(q') + (1 - \alpha) J'_{eq}(q'-1)$$ where $\alpha$ is the weighting factor given on the current gradient $\nabla J_{eq}(q')$. Now, using the average gradient $\nabla J''_{eq}(q')$ in and in place of $\nabla J_{eq}(q')$, we can estimate the SNC factor $\mathbf{G}$ which can then be used to get an estimate of the estimated true image of the $i$-th frame, $\hat{\mathbf{R}}_i$, using . Finally, averaging the estimates for $i=1,2,\cdots,p$ for SNR improvement, the true image is reconstructed as $$\label{speckleless}
\hat{\mathbf{R}} = \sum_{i=1}^p \hat{\mathbf{R}}_i$$

The additive noise in has been ignored so far. However, considering the additive noise in , we get $$\begin{aligned}
\label{noiseAdder}
\hat{\mathbf{R}}_i &= (\mathbf{R}\cdot*\mathbf{U}_i + \mathbf{W}_i)\cdot*\mathbf{G}_i \nonumber\\
                    &= \hat{\mathbf{R}}'_i + \mathbf{W}'_i\end{aligned}$$ where $\mathbf{W}'_i$ is the modified additive noise in the $i$-th frame. From and , we can write $$\label{finalImg}
    \hat{\mathbf{R}} = \hat{\mathbf{R}}' + \mathbf{W}'$$ Due to superior performance of non-subsampled shearlet transform (NSST) to capture the geometric and mathematical properties of an image such as scales, directionality, elongated shapes and oscillations as described in [@shahdoosti2016image], we attempt to denoise the image $\hat{\mathbf{R}}$ with a hard-thresholding on the NSST co-efficients following a similar method as described in [@gupta2014speckle]. Here, $\mathbf{W}'$ is somewhat minimized due to averaging. However, according to [@michailovich2006despeckling], [@achim2001novel], [@garg2018despeckling], the effect of speckle noise is more pronounced compared to the additive noise and hence, the SNR for the additive noise can be considered high in the despeckled image. In this approach, the estimated true image is decomposed into $4$ levels with each having $3$, $3$, $4$ and $4$ directions. The coarse scales of the NSST co-efficients are not thresholded, but the finest scale is hard-thresholded using a tunable low threshold value.

Post-processing
===============

To match with the characteristics of the display monitor and control the overall brightness of the image, further post-processing like gamma correction [@lee2005adjustable] is necessary. The gamma correction [@huang2018speckle] of the estimated true image $\hat{\mathbf{R}}$ is done using $$\label{gamma}
    I(m,n) = \left(\frac{\hat{r}(m,n)}{max(\hat{\mathbf{R}})}\right)^{\gamma}$$ Finally, to control the image contrast, gray level transformation [@suetens2002fundamentals] of the image $I(m,n)$ is done as $$\label{gray}
    G(m,n) =
    \begin{cases}
      0, \hfill \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad  \text{if}\ G(m,n) < W_{low} \\
      \frac{(I(m,n)/max(\mathbf{I})-W_{low})}{(W_{high} - W_{low})}, \\
        \hfill \quad \quad \quad \quad \quad \quad \quad \text{if}\ W_{low} \leq G(m,n) < W_{high}\\
      1, \hfill \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{otherwise}
    \end{cases}$$ where $W_{low}$ and $W_{high}$ are tunable parameters (intended for tuning the contrast of the image) such that $$W_{low} < W_{high} < 1 \nonumber$$ The complete framework of the proposed ultrasound image enhancement and reconstruction process is depicted in Fig. \[Block\].

Results
=======

In this section, the efficacy of our proposed framework for high-resolution B-mode image generation is evaluated on both the simulation and *in-vivo* data. The contents of the paper cover: deconvolution, despecking via MSNE, and post-processing for a complete B-mode image generation. However, as the main contribution is the despeckling algorithm, the simulation study is designed to show the effectiveness of our proposed method for speckle noise estimation and despeckling only. The convergence of the algorithm in the presence of additive noise is also shown to justify the use of the proposed constraint. On the other hand, the *in-vivo* images suffer from low resolution and speckle noise arising from the physical phenomena related with the US image acquisition system. Therefore, according to the signal generation models, the *in-vivo* study includes the $2$-D deconvolution to enhance the resolution as the first phase for all the methods involved for comparing despeckling performance. Two types of deconvolution methods, namely– *b*MCFLMS and cepstrum [@taxt2001two] are investigated. Finaly, post processing stage is included to make a complete investigation of the high-resolution B-mode image generation pipeline in a single research paper.

The quality of the despeckled image is compared with those of SRAD (speckle reducing anisotropic diffusion filter) [@yu2002speckle] and OBNLM (optimized Bayesian non-local means-based filtering) [@coupe2009nonlocal] methods. The performance matrices, used in this case, are SNR (signal-to-noise ratio), PSNR (peak signal-to-noise ratio), SSIM (structural similarity index measure) [@shahdoosti2017two], and NIQE (natural image quality evaluator) [@mittal2013making]. Among these indices SNR, PSNR, SSIM require reference image and hence, cannot be used for *in-vivo* data. In the absence of a reference image, NIQE and visual evaluation are the only ways to evaluate the performance of the proposed algorithm. On the other hand, in case of simulation data with the original noiseless image at hand, we attempt to build the intuition behind different aspects, i.e., the number of image frames to be chosen, runtime, efficacy in preserving small details, of the proposed despeckling algorithm.

Simulation Data
---------------

Simuation data were generated using the ‘Modified Shepp-Logan’ phantom available in MATLAB with size $256\times256$. In the simulation study, we have investigated the despeckling efficacy of the proposed framework and have not considered the PSF effect on the image. Hence, the image was corrupted with speckle noise only as described in [@coupe2009nonlocal] following where $\mathcal{V}(m,n) \sim \mathcal{N}(0,~\sigma^2)$ and as in [@loupas1989adaptive], $\eta = 0.5$ was used. The level of noise was varied by setting $\sigma = \{0.2; 0.4; 0.8\}$. At a particular noise level, the speckle pattern was varied using different $\nu(m,n)$ patterns with the same distribution for different frames. Among the frames used for despeckling using our proposed method, the last frame was despeckled using SRAD and OBNLM for comparison with the proposed algorithm. Here, the implementation platform used were: CPU: Intel$\textsuperscript{\textregistered}$ Core$\textsuperscript{TM}$ i$7$-$8700$K, RAM: $32$ GB, software: MATLAB$\textsuperscript{\textregistered}$, The MathWorks, Natick, MA. Comparing and in additive noiseless case, we can write the speckle pattern of the $k$-th frame as $$\label{speckleNoise}
u_k(m,n) = 1 + r^{-0.5}(m,n)\nu_k(m,n)$$ Therefore, can be used to calculate the true speckle pattern for the simulation data, and subsequently, it can be used to quantify the performance of the proposed despeckling algorithm. As claimed in Section \[sssec:ref1\], increasing the number of image frames in the proposed speckle pattern estimation algorithm has an impact on the accuracy of estimation. In what follows, we attempt to establish a suitable frame number that optimally meets all the assumptions made in Section \[sssec:ref1\], consumes less runtime, and produces a visually pleasant despeckled image. The performance index used in this case is NPM (normalized projection misalignment) defined as $$\begin{aligned}
\mbox{NPM}(q) &= 20\mbox{log}_{10}\left(\frac{\|\rho(q)\|}{\|\mathbf{U}\|}\right)~\mbox{dB} \\
    \rho(q) &= \mathbf{U}-\frac{\mathbf{U}^T\mathbf{\hat{U}}(q)}{\mathbf{\hat{U}}^T(q)\mathbf{\hat{U}}(q)}\mathbf{\hat{U}}(q)\end{aligned}$$ A lower value of NPM indicates better estimation of $\mathbf{U}$. From Table \[tab:NPM\], it is apparent that increasing the level of noise deteriorates the estimation accuracy of the speckle pattern. However, it can be improved by around $6$ dB for the three different noise levels as mentioned above by increasing the number of image frames from $5$ to $20$ in the proposed algorithm. As described in [@xu1995least], for an accurate estimation using the blind multichannel algorithm, the channels should not have common zeros. As we introduce more image frames in the estimation process, the probability of having common zeros decreases. This in turn improves the speckle estimation accuracy. However, increasing the number of image frames imply more computational complexity leading to higher runtime, and at the same time, it causes the violation of the quasi-stationarity assumption for the true image as described earlier in Section \[sssec:ref1\]. The despeckled images using different number of frames for noise level $\sigma = 0.4$ are shown in Fig. \[frames\]. It is apparent from this figure that consideration of frame number greater than $10$ results in visually imperceptible change in the despeckled images. Therefore, in Fig. \[phantom\], we have used $10$ image frames to compare our simulation phantom results with other algorithms.

[|c|c|c|c|c|]{} & &\

------------------------------------------------------------------------

Number of frames & & Runtime (sec)\

------------------------------------------------------------------------

& &\

------------------------------------------------------------------------

$5$ & $-31.65$ & $-25.69$ & $-19.96$ & $0.61$\

------------------------------------------------------------------------

$10$ & $-34.69$ & $-28.72$ & $-22.96$ & $2.01$\

------------------------------------------------------------------------

$15$ & $-36.41$ & $-30.45$ & $-24.68$ & $3.95$\

------------------------------------------------------------------------

$20$ & $-37.64$ & $-31.64$ & $-25.83$ & $6.63$\

\[tab:NPM\]

\[tab:simComp\]

The performnce of the proposed algorithm in comparison to others is illustrated in Fig. \[phantom\]. From Fig. \[phantom\](d), we can deduce that SRAD distorts the texture in the homogeneous region and blurs the small details in the phantom. Again, as evident from Fig. \[phantom\](e), although OBNLM is superior in performance compared to SRAD in preserving texture and edges, it fails to remove speckle noise completely when the noise level is high, e.g., $\sigma = 0.4$. On the other hand, our proposed algorithm (see Fig. \[phantom\](f)) shows significant visual improvement in terms of maintaing original texture, edges and small details compared to the SRAD and OBNLM approaches, and the despeckled image is visually close to the clean phantom image. Again, the quantitative metrics as presented in Table \[tab:simComp\] also demonstrate that our proposed algorithm performs significantly better in terms of quantitative indices SNR, PSNR, and SSIM at different noise levels compared to SRAD and OBNLM. To show the effect of the zero-lag correlation constraint on the convergence profile of the proposed MSNE algorithm, we depict in Fig. \[phantom\](h)-(i) the NPM curve with $20$ dB SNR. As can be seen in \[phantom\](h), the algorithm misconverges near $-16$ dB, whereas in \[phantom\](i), there is no sign of misconvergence and the constrained MSNE algorithm smoothly convergences to around NMP$= -22.96$ dB implying that the zero-lag correlation constraint is effective in preventing misconvergence.

*In-Vivo* Data
--------------

Performance of the proposed complete framework for US B-mode image generation comprising of deconvolution, despeckling and post-processing, respectively, is evaluated on the *in-vivo* data, collected using a commercial SonixTOUCH Research (Ultrasonix Medical Corporation, Richmond BC, Canada) scanner integrated with a linear array transducer, L14-5/38, operating at $10$ MHz with sampling frequency of 40 MHz. These data were collected from the patients who appeared for medical examination at the Medical Centre of Bangladesh University of Engineering and Technology (BUET), Dhaka, Bangladesh. This study was approved by the Institutional Review Board (IRB), and prior patient consent was taken.

The performance evaluation of the deconvolution step is done subjectively as elaborate performance evaluation of the $1$-D *b*MCFLMS algorithm is already done in our published work [@dey2018ultrasonic]. However, the performance of the complete framework is evaluated using two approaches. First, we keep the deconvolution step fixed for all the despeckling algorithms to be compared with the proposed MADS algorithm and evaluate their comparative performance both from visual and quantitative perspectives. Second, we compare subjectively the image generated using our complete framework with the data acquiring machine B-mode image. All the *in-vivo* images shown here were log compressed and dynamic range was set to $35$ dB as described in [@dey2018ultrasonic] for display purpose.

[|c|c|c|]{}

------------------------------------------------------------------------

Data &\

------------------------------------------------------------------------

& Axial & Lateral\

------------------------------------------------------------------------

RF & $0.0379$ & $0.0501$\

------------------------------------------------------------------------

$1$-D deconvolved & $0.0287$ & $0.0497$\

------------------------------------------------------------------------

$2$-D deconvolved & $0.0279$ & $0.0329$\

\[tab:correlation\]

### 2-D Deconvolution Performance Evaluation on *In-Vivo* Data

The successive stages of deconvolution offer images with improved and finer texture as shown in Figs. \[deconv\](a)-(c) and the zoomed-in views of their marked portion in Figs. \[deconv\](d)-(f), respectively. The speckle pattern in Fig. \[deconv\](d) is blurry and highly auto-correlated in the spatial domain as convolution of point speckle with the US PSF results in the spreading of the point spatially and thereby, reduces the resolution of the image. This large spatial coverage of speckle leads to considerable correlation between speckle noise not only in the same frame but also in the consecutive frames leading to a greater number of common zeros. To justify our claim, we select an axial and a lateral line from two consecutive frames along the marked lines in Figs. \[deconv\](a)-(c) and calculate the energy of the normalized correlation among two axial lines as well as lateral lines of two consecutive frames as presented in Table \[tab:correlation\]. Higher value of correlation energy indicates higher correlation among the two frames along the axial or the lateral direction. From Table \[tab:correlation\], it is evident that the raw RF data of a particular frame is more correlated with the next frame in both the axial and the lateral directions compared to those of the $1$-D and the $2$-D deconvolved data. However, the $1$-D deconvolved data has higher lateral correlation compared to that of the $2$-D deconvolved data with nearly the same axial correlation (see Table \[tab:correlation\]). And in Fig. \[deconv\](e), the speckle pattern becomes fiber-like with greater spatial coverage along the lateral direction than the axial direction leading to higher correlation with the next frame in the lateral direction. Finally, in Fig. \[deconv\](f), the lateral correlation is sufficiently removed, and the speckle pattern becomes randomly distributed and uncorrelated. As the speckles in the final $2$-D deconvolved image occupies lesser space, this has additional advantage of reducing common zeros between the speckle patterns in consecutive frames along with increasing resolution of the image. Hence, the speckle pattern estimation can be done efficiently with $5-10$ number of frames without violating the quasi-stationarity assumption of the true image.

### Despeckling Performance on *in-vivo* Data

[|c|c|c|c|]{}

------------------------------------------------------------------------

& OBNLM & SRAD & Proposed MADS\

------------------------------------------------------------------------

RF data- 1 & $7.38$ & $7.69$ & $5.16$\

------------------------------------------------------------------------

RF data- 2 & $7.61$ & $8.49$ & $5.32$\

\[tab:in\_vivo\]

We have used two *in-vivo* breast ultrasound images, RF image-1 and RF image-2, to validate the performance of our proposed algorithm. Similar type of post-processing with $\gamma = 0.97$ in and $W_{low} = 1e-2, W_{high} = 0.98$ in were set for all the images despeckled with different algorithms for illustration purpose. The performance index used here is NIQE that relies on the deviation from the statistical regularities of distortionless images to rate an image as defined in [@shahdoosti2018maximum], [@mittal2013making]. The lower the value of the NIQE metric, the better the quality of the despeckled image. However, the NIQE index was measured on the despeckled image without post-processing. As shown in Table \[tab:in\_vivo\], our proposed algorithm gives the lowest NIQE score and hence, the best quality image compared to that of SRAD and OBNLM. Figs. \[data\_1\] and \[data\_2\] are provided for subjective evaluation of the proposed algorithm. As shown in Figs. \[data\_1\](b) and \[data\_2\](b), SRAD algorithm succeeds in preserving edges although it blurs the texture and degrades the contrast of the image. On the other hand, according to Figs. \[data\_1\](c) and \[data\_2\](c) OBNLM shows superior performance compared to SRAD in preserving undistorted texture and contrast. However, it fails to remove the speckle noise completely from the image. To visually compare the performance of our proposed framework to that of the commercial US image acquiring machine used in this experiment, Figs. \[data\_1\](d)-(e) and \[data\_2\](d)-(e) are portrayed.

From these figures, observe that tissue texture is more prominent in the images provided by our proposed framework compared to those in the machine B-mode images. To facilitate the observations, significant structures of the images are marked with arrow and circle in Figs. \[data\_1\](d)-(e) and \[data\_2\](d)-(e) which show the machine B-mode images have blurred and distorted the tissue structures. Again, the cyst boundary in Fig. \[data\_2\](d) is sharper and well-defined comapred to that of Fig. \[data\_2\](e). The estimated speckle noise patterns of the $5$-th image frame as shown in Figs. \[data\_1\](f) and \[data\_2\](f) contain tissue structures that justify the relevance of true image dependent modeling [@coupe2009nonlocal] of speckle pattern as shown in .

In spite of offering an elegant solution to the speckle removal problem, the proposed framework has a flaw in its complete pipeline as the deconvolution step is not realtime implementable requiring $76$ minutes in total for a single image of $128$ A-lines with each line having $1040$ samples. However, as an alternate approach, the deconvolution step can be replaced by a time-efficient cepstrum-based deconvolution as described in [@taxt2001two] with a little cost paid in image quality as evident from higher NIQE score of $5.59$ (see Fig. \[cepstrum\] (b)) compared to that of $5.32$ (see Fig. \[cepstrum\](a)) using the *b*MCFLMS algorithm, and this brings down the total runtime to $6.3$ seconds. A graphical processing unit (GPU) based deconvolution technique to be investigatd in future may bring down the overall B-mode image generation framework into real-time.

Discussion
==========

In this study, we present a complete framework of signal processing approaches comprising of deconvolution, despeckling, gamma correction, and gray level transformation to produce a high-resolution B-mode image with superior edge and texture from the raw RF image. The parameters for SRAD and OBNLM algorithms were tuned for the lowest NIQE score of the despeckled image. While deriving SRAD and OBNLM algorithms in [@yu2002speckle] and [@coupe2009nonlocal], respectively, deconvolution of raw RF image to enhance resolution was not addressed. Hence, introducing deconvolution prior to SRAD and OBNLM may have resulted in their poor performance. In the proposed framework, the performance of the despeckling algorithm (MADS) is dependent on the number of consecutive image frames to be considered in the MSNE algorithm. As mentioned earlier, there is a tradeoff between the number of frames that can be used without violating the quasi-stationarity assumption of the true US image and the speckle noise estimation accuracy. In our experiment, we observed that five consecutive image frames are good enough for a visually pleasant B-mode image generation. Again, the Lagrange multipliers– $\beta_1$ in the constraint preventing misconvergence of the MSNE algorithm (see ) and $\beta_2$ in the energy constraint of the iterative despeckling algorithm (see ) remain effective once set at an optimum level for a particular US imaging set-up. To make the framework independent of the display monitor, gamma correction as a post-processing step has been introduced. Again, to offer the user a tunable contrast adjustment, two parameters $W_{low}$ and $W_{high}$ have been used in the gray level transformation step.

In addition to providing a guideline for high-resolution B-mode image generation, the paper introduces a method to extract the speckle pattern inherent in a US image. Despite of being a random process, speckle noise is not devoid of information. Since the statistics of the speckle depends on the microstructure of the tissue parenchyma, it can be useful for differentiating between different tissue compositions or types [@wagner1983statistics], [@sehgal1993quantitative].

Conclusion
==========

This paper has dealt with a complete framework for high-resolution ultrasound image reconstruction from raw RF data. The proposed framework relies on SIMO models for both deconvolution and speckle noise estimation, and MISO model for despeckling. In the first step, to enhance the resolution, a $2$-D deconvolution technique has been introduced as an extension of our previously proposed $1$-D *b*MCFLMS algorithm which is necessary prior to despeckling according to the mathematical model of US imaging. In the next step, a novel multiframe-based adaptive speckle noise estimation (MSNE) algorithm estimates the speckle pattern without any *a priori* information on the statistics of the image or the noise pattern. Using the estimated speckle pattern, an energy constrained iterative algorithm estimates the true US image following a MISO model. As the despeckling procedure is completely based on signal generation model and does not involve any kind of ad-hoc filtering operation as reported in the literature, it has resulted in a high quality tissue texture and edges in the image. Finally, gamma correction and gray level transformation have been done as post-processing to complete the pipeline of high quality B-mode image generation. The results have demonstrated the superiority of our proposed despeckling algorithm compared to SRAD and OBNLM methods. Again, the proposed framework offers B-mode image with superior texture and image details compared to those provided by a commercial ultrasound scanner.

As our proposed framework preserves original image features such as texture, details and edges, it may have a far reaching impact on medical imaging for diagnostic purpose. At the same time, the proposed despeckling algorithm may be efficacious in dealing with the speckle noise problem in other imaging such as synthetic aperture radar (SAR) [@wang2019bayesian] and optical coherence tomography (OCT) [@paul2019speckle].

Acknowledgment {#acknowledgment .unnumbered}
==============

This work has been supported by the Higher Education Quality Enhancement Program (HEQEP) of University Grants Commission (UGC), Bangladesh (CPSF-96/BUET/W-2/2017).

[Jayanta Dey]{} received the B.Sc. and M. Sc. degrees in electrical and electronic engineering from the Bangladesh University of Engineering and Technology (BUET), Dhaka, Bangladesh, in 2017 and 2019, respectively. He is currently pursuing his PhD in the BME dept. of Johns Hopkins University.

He was a participant at IEEE Signal Processing Cup 2016, Shanghai, China and his team achieved the 1st position at that competition. His current research interests include digital signal processing, adaptive filtering, image processing, neuroscience, machine learning, and medical imaging.

[Sharmin R. Ara]{} received the BSc degree in Electrical and Electronic Engineering from Bangladesh University of Engineering and Technology (BUET) and MSc degree in Electrical and Computer Engineering from Southern Illinois University, in 2000 and 2007, respectively. She received her Ph.D. degree in Electrical and Electronic Engineering from BUET in 2016. She is working as an assistant professor in the Department of Electrical and Electronic Engineering at East West University, Bangladesh. Her research interests include ultrasonic imaging and tissue characterization.

[Md. Kamrul Hasan]{} received the B.Sc. and M. Sc. degrees in electrical and electronic engineering from the Bangladesh University of Engineering and Technology (BUET), Dhaka, Bangladesh, in 1989 and 1991, respectively, and the M.Eng. and Ph.D. degrees in information and computer sciences from Chiba University, Chiba, Japan, in 1995 and 1997, respectively.

In 1989, he joined BUET as a Lecturer with the Department of Electrical and Electronic Engineering, where he is currently working as a Professor. He was a Postdoctoral Fellow and a Research Associate at Chiba University and Imperial College, London, U.K., respectively. He was a short-term invited Research Fellow at the University of Tokyo, Tokyo, Japan, and a Professor of International Scholars of Kyung Hee University, Seoul, Korea. He has authored or coauthored more than 100 scientific publications. His current research interests include digital signal processing, adaptive filtering, image processing, quantitative ultrasound, elastography, machine learning, and medical imaging. Dr. Hasan is currently an Associate Editor for the IEEE Access.

[^1]: Jayanta Dey, S. R. Ara, and M. K. Hasan are with the Department of Electrical and Electronic Engineering, Bangladesh University of Engineering and Technology, Dhaka, Bangladesh (corresponding author’s e-mail: khasan@eee.buet.ac.bd).
