1 Introduction

Ultrasound (US) imaging is a diagnostic imaging technique to visualize various subcutaneous body structures such as muscles, vessels, and joints, for detecting or diagnosing the pathological symptoms. There are essentially four different modes of scanning as far as US sonography is concerned. Among them, B-mode sonography is more popular due to its extent of applicability and wide range of acceptance among the medical imaging community for diagnostic purpose. The extensive use of US imaging is justified by its non-invasive and non-ionizing nature of characterizing the tissue. The US scanning system sends acoustic signals to probe the nature of tissues in the human body, and based on the reflected pattern, the images are formed. However, the acoustic signals used for probing the tissues result in constructive and destructive interferences causing high- and low-amplitude deflections in the magnitude of the captured data usually coined as speckles [1]. Speckle by its sheer nature is data correlated and cannot be neglected like the other noise interventions [1]. Speckles are not completely noise components though apparently they seem to be. They carry information about the tissue. However, their presence results in spurious analysis by the medical experts. The low-amplitude signals are sometimes treated as blood vessels though they are formed due to the destructive interference of the waves. The formation of the speckles has inspired many researchers to assume the model of formation as multiplicative.

The statistical filters were introduced to effectively address the multiplicative speckled interference in many previous works [2,3,4,5]. The pioneer works in this field are Lee filer [2], Frost filter [3], and Kuan filter [4]. The extent of averaging in these filters depends on the coefficient of variation of the image. The areas with high variance are smoothed at a lesser magnitude compared to the low-variance regions which are categorically smooth in nature. This property in turn preserves the structures apparently. The minimum mean square filter based on the local variance introduced by Weiner [5] is another similar approach in this direction. Nevertheless, these filter duly neglect the noise distribution in the input and the correlated nature of the noise to a large extent. These issues are addressed in their subsequent modifications.

Anisotropic diffusion models inspired by the Perona–Malik model [6] has apparently changed the outlook of the image restoration framework. The Perona–Malik model is a nonlinear diffusion model whose diffusion coefficient is controlled by a nonlinear function of the gradient term.

The models such as speckle reducing anisotropic diffusion (SRAD) [7], oriented SRAD (OSRAD) [8], and direction-preserving anisotropic diffusion (DPAD) [9] are the notable anisotropic diffusion models introduced for despeckling images. In these filters, the diffusion magnitude is controlled by the coefficient of variation which depends on the local mean and variance of the data. For instance, the instantaneous coefficient of diffusion serves as the edge detection function in [7] and the eigenvectors of the structure tensor are used to control the magnitude of diffusion in [8]. The eigenvector in the principal major direction is smoothed less compared to the principal minor direction in case of a 2D matrix. This eventually results in despeckling with less penalization of the structures. There are certain advancements suggested for these filters as well [10, 11] to improve their performance in terms of image restoration. In [10], the authors propose a speckle reducing model where a speckle reducing edge detector is embedded in the well-known “Geodesic snakes” model. Similarly, in [11] the authors introduce a complex diffusion-driven coefficient to drive the process. Recently, another improved diffusion-based method has been proposed in [12] which is an Anisotropic Diffusion Filter with Memory Based on Speckle Statistics (ADMSS) where they have tried to control the diffusion operation on the basis of structural information in memory. However, these filters too ignore the distribution of the data at large. An optimized Bayesian nonlocal means (OBNLM) [13] method is also introduced for speckle reduction which is based on the similarity of the patches of pixels. Later, to eliminate the drawback of over-smoothing and to enable texture preservation, non-local methods with different similarity measures are also introduced, see [14].

Thresholding schemes are employed in large scale in image restoration applications. In particular, multi-resolution models have captured the attention lately and they are extensively employed in image restoration. For instance, wavelet models are prominent in this regard. The wavelet decomposition approximates a scale space representation of the data [15]. The noise also gets represented in the scale space domain, and an appropriate thresholding scheme should effectively remove noise features. Hard and soft (adaptive) thresholding schemes are proposed in the literature to effectively handle the noise interventions. Apart from wavelet models, other advanced versions such as curvelets have also been proposed for restoration activities in the literature, see [16] for details. These methods convert the original image into a logarithmic domain and make use of the Gaussian distribution assumption of the sub-band coefficients. One of the main drawbacks here is related to the inaccuracy in choosing an appropriate threshold. Later, several researchers proposed methods for modified threshold selection process which includes the nonlinear estimator in [17]. In [18], authors apply a log transform on the input to make the data additive in nature but later assumes the wavelet sub-bands to follow non-Gaussian statistics such as heavy-tailed alpha-stable distribution. Recently, a multi-directional 2-D eigenfilter approach for US denoising is proposed by Nagare et al. in [19]. This method uses Translation Invariant Pyramidal Directional filter Bank (TIPDFB) to decompose the image, and thresholding is applied on all the sub-bands to remove noise. All these multiscale methods use the Gaussian distribution assumption in a log domain instead of incorporating speckle characteristics in the denoising process. The additive Gaussian noise has been explored extensively in the literature [20, 21], and a simplified log transformation theoretically converts a multiplicative model to an additive one. However, such oversimplified models tend to neglect the very basic nature of the noise and its distributional characteristics. Even under a log transformation the noise does not seems to become data-uncorrelated [22].

The variational frameworks are the extensions of penalization theory wherein the optimization functional is defined in terms of the regularization and data fidelity aspects of the data. The total variation (TV) regularization is a well-known variational model for image restoration, see [20] for the details. However, the linear approximation of the model results in undesired effects in the restored output. Moreover, the model assumes a Gaussian distributed white noise intervention in the data. A notable modification to this model to incorporate the multiplicative Gamma noise is introduced in [22]. Here, a Bayesian formulation is explored for reinterpreting the minimization problem as a posterior probability minimization model. The model is found to perform well in case of multiplicative speckles following a Gamma law. Many variant thoughts in this direction are seen in works [23] and [24]. However, these methods too are not efficient enough to maintain small textures and details present in the data, even though they ensure the preservation of edge details. To address this issue to a considerable extend, a new category of variational methods was introduced, which is based on the non-local (TV) framework proposed by Gilboa and Osher in [25]. Some variations of this model are found in [26] and [27], where the later works propose a nonlocal total bounded variation regularization model for speckle removal. Compared to the other methods, the non-local (TV) algorithms have succeeded in ensuring the efficiency in retaining small details and textures. In [28], authors proposed a combination of diffusion and TV filters where a phase asymmetry measure for edge detection is being used to distinguish edges and constant intensity regions.

Many of the models discussed so far seek to find a despeckling solution that eventually retains the edges and other details present in the data to a considerable extend. Nevertheless, the images captured under the US modality are observed to be deficient in contrast aspects as well. An uneven contrast distribution duly hinders the subsequent analysis phase, thus spoiling the diagnosis and detection phase of automated systems. This has been a matter of concern among the research community. The untiring efforts in addressing the uneven distribution of the data resulted in Retinex-based models, which are extensively used in low-contrast photography reconstruction algorithms.

The Retinex theory has evolved form the Land’s key discovery of colour constancy phenomenon in human visual system [29]. This study explained that the human eye can adjust with illumination changes and it can correctly distinguish and match colours of objects even in varying illuminations, so we have

$$\begin{aligned} U(x,y)= L(x,y) R(x,y), \end{aligned}$$
(1)

where \(R\in [0,1],\) \(L\in \mathbb {R}^2,\) and \(U\in \mathbb {R}^2.\) It suggests that any pixel (x,y) in an intensity image U can be represented as a product of reflectance R and illumination L. Readers are invited to refer [29]. This decomposition of images helps in eliminating the effect of nonuniform illumination. However, estimation of illumination from an intensity image is an ill-posed problem.

Kimmel et al. in [30] introduced a variational model for Retinex theory. The authors described this problem as a quadratic optimization equation where \(L_{2}\) norm of the illumination and reflectance part is evaluated to ensure the smooth and visually appealing result. The optimization model is given by:

$$\begin{aligned} \min E[l]=\int _\varOmega \left( \Vert \nabla l\Vert _2^2 +\alpha (l-u)^2+\beta \Vert \nabla (l-u)\Vert _2^2\right) \mathrm{d}x\mathrm{d}y, \end{aligned}$$
(2)

where l and u represent the illumination and intensity matrices in log domain and \(u=l+r\). However, smoothing the reflectance matrix results in penalization of textures and details. To address this drawback, a (TV) Retinex model is introduced in [31].

$$\begin{aligned} \min _{l\in L^2 (\varOmega ) ,r \in L^1(\varOmega )} E[l,r]&=\int _\varOmega \Vert \nabla r\Vert \mathrm{d}x\mathrm{d}y+\alpha \int _\varOmega \Vert \nabla l\Vert _2^2 \mathrm{d}x\mathrm{d}y\nonumber \\&\quad +\beta {\int _\varOmega } (l-r-u)^2 \mathrm{d}x\mathrm{d}y+\mu \int _\varOmega l^2 \mathrm{d}x\mathrm{d}y. \end{aligned}$$
(3)

In this model, the authors use (TV) norm in the reflectance term r for preserving more details, see equation (3). Later, inspired from the non-local (TV) proposed by Gilboa and Osher in [25], a non-local Retinex version was also introduced in [32], see [25] and [32] for the details of non-local operators. All of these variational Retinex methods are proposed for enhancing natural images by reducing non-uniform illumination. In addition to this, a few variants of variational Retinex models are also proposed for image enhancement in some specific domains, see [33].

As studied in many related works, the US images are corrupted by data-correlated multiplicative noise. The noise in many of the instances is observed to follow a Gamma [22] distribution. A model that duly respects the distribution characteristic of the speckle noise is expected to perform better than its counterparts in restoring the data. Moreover, most of the US images used in various image analysis applications are of poor quality due to various factors influencing their formation as discussed earlier in this paper. These facts call for a holistic model which can address all the degradations simultaneously. For instance, the contrast unevenness of the US data duly hinders the analysis phase and sometimes even results in spurious diagnosis. Most of the despeckling procedures in the literature at times do not considerably improve the contrast and resolution aspects of the data well. This is another concern which we intent to address in this study. The attempt to address both denoising and enhancement by a variational Retinex model is relatively a new approach. A similar attempt for Gaussian distributed noise is found in [34] and [35], and along the similar lines, the model is further extended to Poisson noise in [36]. This paper introduces a model that does not neglect noise distribution and restores and enhances the US images without penalizing the details present in them. A non-local variational Retinex framework is used to design the model.

The rest of the article is organized as follows. In Sect. 2, the proposed framework for denoising and enhancement is explained. Numerical implementation of the proposed optimization problem is also elaborated in this section. The details of the experimental analysis conducted are given in Sect. 3. The last section provides concluding remarks and mentions the scope of future work.

2 Proposed framework

As studied in many previous works, a speckle distorted image is generally represented as

$$\begin{aligned} U_0=U\times n, \end{aligned}$$
(4)

where n is multiplicative data-correlated noise, U denotes the original amplitude data (unobserved one), and \(U_0\) represents its speckled version. The noise distribution is observed to follow a Gamma law. This claim has been substantiated in many previous works [22, 23, 37]. Now, we derive the optimization model for denoising using the Bayesian framework as done in many preceding works, see [22]. The Bayesian law says

$$\begin{aligned} P(U|U_0)= \frac{P(U_0|U) P(U)}{P(U_0)}, \end{aligned}$$
(5)

where \(P(U|U_0)\) is the posterior probability (i.e. conditional probability of U given \(U_0\)), \(P(U_0|U)\) is the likelihood probability, and P(U) is the prior probability on U. Maximizing the posterior probability gives the restored version of the image, so we have:

$$\begin{aligned} \max _U P(U|U_0)=\max \left\{ P(U_0|U) P(U)\right\} , \end{aligned}$$
(6)

since \(P(U_0)\) is constant with respect to U,  we omit this term from further evaluations. However, maximizing the function and its log likelihood is one and the same (as log function is monotonic) we have

$$\begin{aligned} \max _U \log (P(U|U_0))=\max \left\{ \log (P(U_0|U))+\log (P(U))\right\} . \end{aligned}$$
(7)

Further maximizing a function is the same as minimizing its negative, so

$$\begin{aligned} \min _U \{-\log (P(U|U_0))\}=\min _U \left\{ -\log (P(U_0|U))-\log (P(U))\right\} . \end{aligned}$$
(8)

According to our assumption, the noise n follows a Gamma law and based on the proof given in [22], it implies

$$\begin{aligned} P(U_0|U)=\frac{J^J}{U^{J}\varGamma (J)}U_0^{J-1}e^{-\frac{JU_0}{U}}, \end{aligned}$$
(9)

where \(\varGamma \) is the Gamma function and J represents the number of images used in averaging; it is one in this case as we are using a single image. We also assume U to follow a Gibbs prior, i.e. \(P(U)=e^{-\lambda \phi (U)},\) where \(\phi (.)\) is the weighted and weberized non-local TBV prior. Therefore, equation (8) can be represented as

$$\begin{aligned} \min _U \{-\log (P(U|U_0))\}=\min _U \left\{ \log (U)+\frac{U_0}{U}+\lambda _x \phi (U)\right\} , \end{aligned}$$
(10)

where \(\lambda _x\in (0,\infty )\) denotes the positive scalar regularization parameter. Considering \(\varOmega \) as a image domain, the derived functional term takes the form:

$$\begin{aligned} F(U)=\int _\varOmega \log (U)+(U_{0}/U)+ \lambda _x \int _\varOmega \phi (U) \mathrm{d}x\mathrm{d}y. \end{aligned}$$
(11)

Unlike the model in [22], we use a novel weighted and weberized non-local TBV prior to obtain a better denoising and contrast improvement to the input images. Based on these assumptions, the model proposed herein is represented as:

$$\begin{aligned} \min _{l,r,u}\{E(l,r,u)\}&=\lambda _0 \int _\varOmega \Vert \nabla l\Vert _2^2 \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _1 \int _\varOmega W(r) \Vert \nabla \hat{r}\Vert \mathrm{d}x\mathrm{d}y\nonumber \\&\quad +\lambda _2 \int _\varOmega (exp(r)-1/2)^2 \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _3 \int _\varOmega (W(u) \Vert \nabla \hat{u}\Vert +\beta \Vert u\Vert _2^2) \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _4 \int _\varOmega (r-u+l)^2 \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _5 \int _\varOmega \log (U)+(U_{0}/U), \end{aligned}$$
(12)

Inspired by previous variational methods, this model also evaluates the \(L_{2}\) norm of illumination \((|\nabla l\Vert _2^2)\) to ensure its smoothness. As the reflectance part of the image is containing more texture information, it is proposed to use a weighted non-local (TV) \((W(r) \Vert \nabla \hat{r}\Vert )\). As mentioned by Kimmel et al., in [30], the reflectance is in the range [0, 1] and it is forced to be close to average value by adding a term \((exp(r)-1/2)^2\) in minimization function. The term \((r-u+l)\) ensures that the intensity matrix evaluated in each step is a combination of corresponding r and l estimates and vice versa. In this proposed model, we incorporate a weighted non-local total bounded variation (TBV) [27] of intensity matrix \((W(u) \Vert \nabla \hat{u}\Vert +\beta \Vert u\Vert _2^2)\) to eliminate the noise in images. The TBV term used is a combination of NLTV and \(L^2\) norm of the intensity. As the value of \(\beta \) increases, the smoothing of u also increases. Hence, TBV is efficient in giving better results when the data are heavily corrupted by noise. The drawback of using \(L^2\) norm is that it does not care about the edges and details in an image, which is duly addressed by NLTV term in TBV. The intensity smoothing is constrained by the term \(\int _\varOmega \log (U)+(U_{0}/U)\) which denotes the data fidelity. The term is conditionally convex (i.e. convex if \(2U_0>U.\)), and it is derived using the Bayesian MAP estimator under the assumption of a Gamma distributed noise in the input image, refer equation (10) and (11). This term tends to retain more textures and small details while reducing the speckle and eventually enhancing the input. The reflectance and intensity are assumed to be the properties of the object being imaged, and hence, they are expected to carry non-smooth information. Therefore, the non-local (TV) norm is employed to ensure proper retention of the structures. Furthermore, \(\hat{u}\) and \(\hat{r}\) denote the Weberized intensity and reflectance terms, used to improve the visual quality of the resultant image, see [38] for the details. To further assure the preservation of fine details, use of weight functions W(u) and W(r) is proposed in this model. The weight W(u) is evaluated as \(1+\frac{1}{1+\Vert \nabla u\Vert _{1}},\) where weight W(u) increases as the magnitude of gradient decreases and reaches a maximum value two, for a homogeneous region. Similarly, W(r) is also estimated. The luminance being the property of the source, an \(L^2\) norm is employed to ensure its smooth transition. The control parameters \(\lambda _0\) to \(\lambda _5\) can take positive scalar values which can be tuned to get the optimally enhanced and denoised result. Since the optimization framework contains six different terms, the solution is not trivial and it is hard to employ a simultaneous solution for all the terms. As one term influences the other one, we need to seek for a solution by separating the variables, so a split scheme is usually employed for the same.

The optimization problem in equation (12) shall be split into three separate minimization problems (in terms of the three quantities r,  l and u) as given below:

$$\begin{aligned} \min _{r}\{E(r)\}&= \lambda _1 \int _\varOmega W(r) \Vert \nabla \hat{r}\Vert \mathrm{d}x\mathrm{d}y \nonumber \\&\quad +\lambda _2 \int _\varOmega (exp(r)-1/2)^2 \mathrm{d}x\mathrm{d}y\nonumber \\&\quad + \lambda _4 \int _\varOmega (r-u+l)^2 \mathrm{d}x\mathrm{d}y , \end{aligned}$$
(13)
$$\begin{aligned} \min _{u}\{E(u)\}&= \lambda _3 \int _\varOmega (W(u) \Vert \nabla \hat{u}\Vert +\beta \Vert u\Vert _2^2) \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _4 \int _\varOmega (r-u+l)^2 \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _5 \int _\varOmega \log (U)+(U_{0}/U) \mathrm{d}x\mathrm{d}y, \end{aligned}$$
(14)

and

$$\begin{aligned} \min _{l}\{E(l)\}&=\lambda _0 \int _\varOmega \Vert \nabla l\Vert _2^2 \mathrm{d}x\mathrm{d}y \nonumber \\&\quad + \lambda _4 \int _\varOmega (r-u+l)^2 \mathrm{d}x\mathrm{d}y. \end{aligned}$$
(15)

These equations are to be solved iteratively to get the final result.

2.1 Numerical implementation

There are various alternating minimization schemes for solving the above-mentioned problem. The numerical implementation seeks for a solution which converges at a higher rate. Given many such implementations such as alternating method [39], projection method [40], primal-dual method [6], and Bregman method [41], we chose the split-Bregman solution for the reason that it provides a stable solution with a high convergence rate and less parameter sensitivity. Therefore, the above equations (13), (14), and (15) are solved using the split-Bregman numerical optimization technique, see [41] and [42] for more details. With regard to the split-Bregman iteration (see [41]), we introduce new constraints as \(a=\nabla \hat{r}\) and \(b=\nabla \hat{u}\) along with the auxiliary variables \(d_{1}\) and \(d_{2}\), and then, equations (13) and (14) get transformed as

$$\begin{aligned} \min _{r}\{E(r)\}&= \lambda _1 W(r) \Vert a\Vert \nonumber \\&\quad +\lambda _2 (exp(r)-1/2)^2 \nonumber \\&\quad + \lambda _4 (r-u+l)^2\nonumber \\&\quad + \lambda \Vert a-\nabla \hat{r}-d_{1}\Vert ^{2}_{2} \end{aligned}$$
(16)

and

$$\begin{aligned} \min _{u}\{E(u)\}&= \lambda _3 (W(u) \Vert b\Vert + \beta \Vert u\Vert _2^2) \nonumber \\&\quad + \lambda _4 (r-u+l)^2 \nonumber \\&\quad + \lambda _5 \log (U)+(U_{0}/U)\nonumber \\&\quad +\alpha \Vert b-\nabla \hat{u}-d_{2}\Vert ^{2}_{2}. \end{aligned}$$
(17)

The above equations can be further split as given below (we assume that \(\lambda _1=1\) and \(\lambda _3=1\)):

$$\begin{aligned} r^{k+1}&= \min _{r}\Big \{\lambda _2 (exp(r)-1/2)^2 + \lambda _4 (r-u+l)^2\nonumber \\&\quad + \lambda \Vert a-\nabla \hat{r}-d_{1}^{k}\Vert ^{2}_{2}\Big \}, \end{aligned}$$
(18)
$$\begin{aligned} a^{k+1}&=\min _{r}\Big \{ W(r) \Vert a\Vert +\lambda \Vert a-\nabla \hat{r}-d_{1}^{k}\Vert ^{2}_{2}\Big \}, \end{aligned}$$
(19)
$$\begin{aligned} u^{k+1}&= \min _{u}\Big \{\lambda _5 (\log \ U+(U_{0}/U)) + \beta \Vert u\Vert _2^2+ \lambda _4 (r-u+l)^2 \nonumber \\&\quad + \alpha \Vert b-\nabla \hat{u}-d_{2}^{k}\Vert ^{2}_{2}\Big \}, \end{aligned}$$
(20)

and

$$\begin{aligned} b^{k+1}&=\min _{u}\Big \{ W(u)\ \Vert b\Vert +\alpha \Vert b-\nabla \hat{u}-d_{2}^{k}\Vert ^{2}_{2}\Big \}, \end{aligned}$$
(21)

The auxiliary variables \(d_{1}\) and \(d_{2}\) get updated after each corresponding Bregman iterations as follows:

$$\begin{aligned} d_{1}^{k+1}&= d_{1}^{k}+(\nabla \hat{r}-a^{k+1}). \end{aligned}$$
(22)

and

$$\begin{aligned} d_{2}^{k+1}&= d_{2}^{k}+(\nabla \hat{u}-b^{k+1}). \end{aligned}$$
(23)

The Euler–Lagrange is applied on equation (18) and (20) to evaluate the derivative of the functional, which is later solved using the Fourier transform. The iterative procedure for the restoration is given in Algorithm 1:

Fig. 1
figure 1

Restoration results of US B-mode input a noisy input, b SRAD, c DPAD, d OBNLM, e NLTBV, f ADMSS, g PFDTV h eigenfilter, i proposed method

Full size image
Fig. 2
figure 2

Restoration results of US B-mode input a noisy input, b SRAD, c DPAD, d OBNLM, e NLTBV, f ADMSS, g PFDTV h eigenfilter, i proposed method

Full size image
Fig. 3
figure 3

Restoration results of US B-mode input a noisy input, b SRAD, c DPAD, d OBNLM, e NLTBV, f ADMSS, g PFDTV, h eigenfilter, i proposed method

Full size image
Fig. 4
figure 4

Restoration results of US B-mode input a noisy input, b SRAD, c DPAD, d OBNLM, e NLTBV, f ADMSS, g PFDTV, h eigenfilter, i proposed method

Full size image
Fig. 5
figure 5

Restoration results of US B-mode input a noisy input, b SRAD, c DPAD, d OBNLM, e NLTBV, f ADMSS, g PFDTV, h eigenfilter, i proposed method

Full size image
Fig. 6
figure 6

Restoration of synthetic noisy data (a) input image after adding Gamma noise of variance=0.25; (b) SRAD; (c)DPAD; (d) OBNLM; (e) NLTBV; (f) ADMSS; (g) PFDTV; (h) eigenfilter; (i) proposed method

Full size image
Fig. 7
figure 7

Restoration of a 1D signal, sub-figures (b)-(h) contains original signal (blue), input signal after adding gamma noise of variance=0.25 (red) and different restored signals (green): (a) original signal, (b) SRAD, (c) DPAD, (d) OBNLM, (e) NLTBV, (f) ADMSS, (g) PFDTV, (h) eigenfilter, (i) proposed method

Full size image
Fig. 8
figure 8

Parameter analysis using different quality measures (QM) Row1: change in QM with respect to \(\lambda \); Row2: change in QM with respect to \(\alpha \); Row3: change in QM with respect to \(\lambda _2\); Row4: change in QM with respect to \(\lambda _5\); Row5: change in QM with respect to \(\lambda _4\); Row6: change in QM with respect to \(\beta \); Row7: change in QM with respect to \(\lambda _0\)

Full size image
Table 1 Comparison on the basis of different quality metrics
Full size table
figure c

3 Experimental results

For the experimental evaluation, we use the US image dataset provided by Signal Processing Laboratory (http://splab.cz/en/download/) and COVID-19 US data provided in [43]. The former dataset includes B-mode US images of the common carotid artery (CCA). The experimental studies are conducted using MATLAB R2018b on a system having Intel(R) Core(TM) i7-9750H CPU @2.60GHz processor with 16GB of RAM. For comparative study, we use the popular despeckling algorithms such as SRAD [7], DPAD [9], OBNLM [13], NLTBV [27], ADMSS [12], PFDTV [28], and eigenfilter [19].

The restoration results obtained for various comparative methods are given in Figs. 1, 2, 3, 4, and 5 for a visual comparison. In Fig. 1, the details present on the upper portion of the image are missing in most of the restoration results except the proposed method. The performance shown by the eigenfilter is comparatively better than the others; however, the enhancement capability of the proposed algorithm supersedes the others. Further, the proposed model is observed to retain small details present in the image pretty well. A similar kind of performance is visible from Figs. 2 and 3 as well, where the details on the upper portion of the image are least preserved by the other methods including OBNLM. The restoration results obtained using eigenfilter ([19]) appears low in contrast and blurry in appearance, whereas the proposed algorithm gives a better contrast enhancement and denoising at the same time, which makes the model superior to the other comparative models. In Fig. 4, most of the comparing methods except the eigenfilter did not succeed in maintaining textures. Moreover, the details present in the dark regions are entirely washed out in all the comparative results. In contrast, the proposed method has promisingly enhanced the dark regions and preserved details from being ignored in further steps. The textures present in tissues are also preserved in the line of denoising. A similar performance is visible in Fig. 5 as well.

Performance comparison of a synthetic image is depicted in Fig. 6. The noisy input is given in sub-figure (a). As observable from the results, the presence of speckle noise has degraded the contrast of the image as well. Most of the diffusion-based restoring algorithms such as SRAD, DPAD, and ADMSS give poor contrast enhancement, see sub-figures (b), (c), and (f), respectively. The OBNLM has considerably done well in maintaining the contrast, but as shown in Fig. 6d, along with the smoothing (due to despeckling) some artefacts are also introduced in the restored images. Similarly, NLTBV method also restores the original contrast of the image to some extent. The denoising achieved by this model in low-intensity areas is remarkable; however, in high-intensity areas, the performance is moderately degraded. The restoration provided by eigenfilter method is good in suppressing the noise and restoring the contrast aspects of the data, but the resultant image appears blurred with minimal sharp details or textures. On the contrary, the PFDTV result gives sharp edges and details, but contrast between the dark pixel regions like the back square and the background is very low. Comparatively, the method proposed herein gives a better performance even in case of highly corrupted data by giving a consistent performance in low- and high-intensity regions, see sub-figure (h). The method preserves edges and other sharp details efficiently, and the contrast improvement achieved eventually improves the analysis of real scan-images, as it reveals the hidden details which are otherwise unobserved. The same performance pattern is visible in restoring a one-dimensional signal too, refer Fig. 7. Here, sub-figures (b), (c), and (f) represent the SRAD, DPAD, and ADMSS results, respectively, where the restored signal (Green) appears noise-free, but the amplitude or pixel values are getting diminished in the restored data. The range of output amplitude is also reduced to [0-50], which indicates a degradation in contrast. Moreover, unlike the noisy input signal (Red), the shifts from different pixel levels are not sharp in the output, and the lines are mostly straight, which indicates the poorly preserved edges and textures. The OBNLM result given in sub-figure (d) is not shown to degrade the pixel values, but the denoising performance is considerably low. The sub-figure (e) provides the result of NLTBV, where the denoising of high pixel values (corresponds to bright areas in a 2D image) is compromised. Nevertheless, denoising of the rest of the parts appears comparatively better. The sub-figure (g) shows the PFDTV result, where it succeeds in preserving sharp edges, but the black region is not distinguishable from the background. The eigenfilter result given in sub-figure (h) succeeds in removing noise; however, texture preservation capability seems low. The output pixel values are slightly reduced in this result, especially in high-intensity regions. The proposed result is given in sub-figure (h), which improves the contrast of the output by enhancing the intermediate intensity levels, which in turn makes the darker regions bright for better visibility. As speckle noise contains useful information, the textures in noisy input (red) are well preserved in this result and edges are also properly retained.

Various quality metrics are used in this study to analyse the performance of the proposed method with different comparing algorithms. These metrics are chosen to evaluate the efficiency of the system in terms of various parameters, such as noise reduction and contrast enhancement. The equivalent number of looks (ENL), naturalness image quality evaluator (NIQE), entropy, global contrast factor (GCF), and contrast-to-noise ratio (CNR) are the different quality measures used in this study. The parameter values used for the experiments are \(\lambda =1.4, \alpha =0.2, \lambda _0=0.1, \lambda _2=8, \lambda _4=800, \lambda _5=0.001, \beta =0.2\). These parameters are set based on various experimental studies performed with different image data. An adaptive selection of them is still an open problem that is being investigated. The detailed analysis is conducted based on these quality metrics to fix the parameter value which is given in Fig. 8. Different rows in this figure represent the influence of each parameter on the quality metrics under study. In the first row, the effect of parameter \(\lambda \) is given. The quality metrics ENL, CNR, and NIQE increase with the increase in \(\lambda \), whereas entropy and GCF decrease. Hence, to obtain the optimal performance according to all quality measures, we have chosen \(\lambda =1.4\). Similarly, all other parameters are also tuned based on the corresponding rows in Fig. 8.

The ENL is a well-known blind quality metric (the measure which does not demand the original data for its computation) used to evaluate the performance of despeckling filters [44]. The higher value of ENL indicates a better speckle removing capability. We have evaluated the same for different region of interests (ROI). The ENL value obtained for one particular ROI is given in Table 1. The methods such as eigenfilter, PFDTV, and DPAD give comparatively good ENL values than the other diffusion and variational methods used for the comparison. However, the ENL metric evaluated for the proposed method is considerably higher for all input images. The NIQE [45] is another widely used blind quality evaluator. It is based on a natural scene statistical model, a low value of which indicates less distortion due to noise or any other sources. In general, NIQE is efficient in detecting distortions present in an image. As per the tabulated NIQE values in Table 1, we notice that the comparing methods like OBNLM, DPAD, and NLTBV are giving less distorted results than other recent methods. Nevertheless, the proposed method is giving a comparatively lower NIQE values in all the test cases, which indicates restoring ability of the model.

Other than denoising capability, the detail-preserving ability is also evaluated using the metric: entropy [46]. A high entropy indicates an increased amount of details present in an image. The comparing methods like OBNLM, NLTBV, and eigenfilter are giving a slightly higher performance in detail preservation. In case of the proposed method, the tabulated entropy values are higher than the comparing ones for all the test images, which implies that the illumination correction along with restoration is useful to preserve small hidden details. We have also used the GCF [47] for studying the contrast improvement. As we can infer from the Table 1, the contrast measure obtained by OBNLM is close to the proposed method in some test cases, but the details present in this image are categorically compromised. However, the proposed method shows a consistent global contrast improvement in all the test cases.

The CNR [48] is similar to signal-to-noise ratio (SNR), and it is used to evaluate the quality of restoration. A high CNR value indicates better denoising and a contrast enhancement. We have evaluated the same on different ROIs, and one such result is included in Table 1. This tabulated value also supports the superior restoration capability of the proposed method. Moreover, the average execution time taken by this method for \(500 \times 500\) image is 6.537 s, which supports the use of this algorithm on real-time data.

4 Conclusion

A weberized non-local TBV Retinex-based restoration algorithm for despeckling has been proposed in this paper, which does the denoising along with contrast enhancement and illumination correction. The use of weighted and weberized non-local TV and TBV norms aids in retaining fine textures and improving the visual contrast. The visual and quantitative evaluation conducted here votes for the better denoising and detail-preserving capacity of the proposed algorithm in addition to its contrast enhancement ability. The time taken for the execution is in seconds, which makes it suitable for real-time applications. The adaptive parameter selection for getting the optimal result is an open problem, and this is a problem to be investigated in forthcoming research paper.