Keywords

1 Introduction

Images are always degraded by noises in process of formulation, transmission and preservation. And in many applications, there is more than one kind of noises. For example, in the low-photon-counting systems like micro focus X-ray detection, astronomy and fluorescence microscopy, only limited photons are collected because of system requirement or physical constraints. Usually, the noise is modelled as Poisson distribution. And the electronic fluctuations and intrinsic thermal are always assumed as additive Gaussian noise [1]. Thus, Poisson noise and Gaussian noise co-exist in such systems. It’s essential to remove the mixed noise for better application. There are some but not many methods for mixed Gaussian-Poisson noise removal. PURE-LET [1], GAT+BM3D [2], GAT+BLS-GSM [3] are considered state-of-the-art. They are almost based on transform domain like wavelet.

The most popular noisy image reconstruction method should be the compressive sensing based methods. In the compressive sensing framework, the original clean image is sparse while the noisy one is not. Thus, the denoising process is actually to recover the sparsity of the image [4], i.e. to decompose the image into a linear combination of a small number of dictionary atoms [5]. It is also called synthesis sparsity. The synthesis sparsity explores the local information of the images and has shown great results. There are two main problems in sparse reconstruction. One is dictionary learning, the other is sparse decomposition [6]. The most representative sparse method should be K-SVD [7]. To improve the sparse reconstruction quality, it is gaining growing attractions to combine it with more image priors. Nonlocal similarity is the most common indicating that there are many similar patches in the images even if they are not in a local and adjacent area. Researches show that nonlocal similarity is good for reconstructing image structures [2, 8]. Actually, clustering-based dictionary learning like K-LLD [9] is an indirect way to use the non local similarity in sparse reconstruction. And in recent years, the CSR [10] proposed by Dong et al. has been widely considered, which combines nonlocal similarity and sparsity to get the double sparse constraints that help to preserve more details in reconstruction. Another literature [11] combines ideas from nonlocal sparse models and the Gaussian scale mixture model to keep the sharpness and suppress undesirable artifacts. However, the methods mentioned above are aimed at single noise model. And also, sparsity itself is limited because images can only lie in a low dimensional subspace. Given this, a dual analysis viewpoint to sparse representations called cosparsity or analysis sparsity was proposed and has attracted more and more attentions. Cosparsity shifts its focus from the nonzero decomposition coefficients of synthesis dictionary to the zeros in analysis dictionary [12]. The differences between these two models were compared in detail in [12]. Nevertheless, unlike synthespis sparsity, the understanding of analysis sparsity remains shallow and scarce today. Among them, [13] introduced many greedy-like methods in cosparse frame-work including Greedy Analysis Pursuit (GAP), analysis CoSaMP (ACoSaMP), Analysis SP (ASP), etc. And Rubinstein et al. proposed the dual K-SVD based on cosparsity called AK-SVD [14] in 2013, which outperformed K-SVD in removing Gaussian noise. In [15] the cosparse model is cast into row-to-row optimizations and use quadratic programming to obtain sparseness maximization results. As the sparsity and cosparsity are complementary, it’s reasonable to combine them together to promote reconstruction quality. In this paper, we propose a new reconstruction method for images degraded by mixed Gaussian-Poisson noise based on sparsity, cosparsity and nonlocal similarity. And in the follow-up to this article, we refer ‘sparse/sparsity’ to synthesis case and ‘cosparse/cosparsity’ to analysis one.

2 The Mixed Gaussian-Poisson Noise Model

Assume that \(y\in \mathbb {R}\) is corrupted by mixed Gaussian-Poisson noise, then the observation model is as follow:

$$\begin{aligned} y=Poisson(u)+n, \end{aligned}$$
(1)

where Poisson(u) means that the original image \(u\in \mathbb {R}^N\) is corrupted by Poisson noise; \(n\sim \mathcal {N}(0,\sigma ^2I)\) is the additive Gaussian noise; N is the product of image length and width.

Usually, the probability density function (PDF) based on joint probability distribution for formula (1) is complicated and thus the corresponding objective function is difficult to solve. One of the common ways to tackle it is to transform the mixed noise into a Gaussian one using Generalized Anscombe Transform (GAT) [16]. Then the problem turns into the additive case that is easier to solve. While in paper [17], the authors built up the PDF based on independent probability distribution to simplify the objective function in a dual adaptive regularization (DAR) scheme to remove the mixed noise. In this paper, we adopt the same strategy as [17]. Under the independent probability distribution, the PDF changes into:

$$\begin{aligned} P_{mixed}(y|u)=\prod _{k=1}^N\dfrac{(u_k)^{y_k}e^{-u_k}}{y_k!}\dfrac{1}{\sqrt{2\pi \sigma ^2}}e^{-\dfrac{\left\| y-u\right\| _2^2}{2\sigma ^2}}, \end{aligned}$$
(2)

where \(u_k\) and \(y_k\) is the \(k_{th}\) component of u and y.

Then the log-likelihood of formula (2) under MAP criterion is:

$$\begin{aligned} \begin{aligned}T_0(u)&=\min \limits _{u}-\sum \limits _{k=1}^N\log \dfrac{(u_k)^{y_k}e^{-u_k}}{y_k!}-\log \dfrac{1}{\sqrt{2\pi \sigma ^2}}e^{-\dfrac{\left\| y-u\right\| _2^2}{2\sigma ^2}}\\&=\min \limits _{u}\sum \limits _{k=1}^N(u_k-{y_k}\log {u_k})+\left\| y-u\right\| _2^2.\end{aligned} \end{aligned}$$
(3)

Motivated by [18], the Taylor’s approximation is used to transform the logarithmic term into a quadratic term. Therefore, the objective function becomes: at \(i_{th}\) iteration:

$$\begin{aligned} T_0(u)=\min \limits _{u}\left\| u-s^i\right\| _2^2+\frac{2}{\eta _i}\left\| y-u\right\| _2^2, \end{aligned}$$
(4)
$$\begin{aligned} s^i=u^i-\frac{1}{\eta ^i}\nabla Poi(u^i), \end{aligned}$$
(5)

where \(\eta ^i\in \mathbb {R}\) is the second-order coefficient and is calculated via Barzilai-Borwein [19]; \(\nabla \) is the gradient operator; \(Poi(u)=\sum \limits _{k=1}^N(u_k-{y_k}\log {u_k})\) is the Poisson component.

3 The Proposed Method

We will formulate our objective function in this section on the basis of Sect. 2. Firstly, we use a sparse melioration based on nonlocal similarity to further constrain the synthesis sparse coefficients. Then, the mixed Gaussian-Poisson data fidelity term and cosparse constraint are combined with the modified synthesis model.

$$\begin{aligned} T(u)=\min \limits _{\alpha _m,u}\sum \limits _{m=1}\left\| R_mu-D\alpha _m\right\| _2^2+\mu \left\| \alpha _m-\beta _m\right\| _1 \nonumber \\ \ \ \ \ \ \ \ \ +\,\tau (\left\| u-s^i\right\| _2^2+\frac{2}{\eta _i}\left\| y-u\right\| _2^2)+\lambda \left\| \varOmega u\right\| _0. \end{aligned}$$
(6)

Thus, the objective function takes use of synthesis sparsity, analysis sparsity and nonlocal similarity and combines these priors together to improve the reconstruction quality.

To make it more intuitive, we rewrite the objective function (6) as:

$$\begin{aligned} T(u)=\min \limits _{\alpha ,u}\left\| u-D\alpha \right\| _2^2\,+\,\mu \left\| \alpha -\beta \right\| _1\,+\,\tau (\left\| u-s^i\right\| _2^2\,+\,\frac{2}{\eta _i}\left\| y-u\right\| _2^2)\,+\,\lambda \left\| \varOmega u\right\| _0, \end{aligned}$$
(7)

where \(R_m\) in (6) is the operator to extract a patch in \(\mathbb {R}^n\) at location m; \(\eta ,\tau \) and \(\lambda \) are parameters to balance the data fidelity term and constraint terms; \(D=[D_1;D_2;\cdots D_c]\in \mathbb {R}^{n*c}\) is the synthesis dictionary. We used the same training strategy in [20] to get a clustering dictionary. \(\varOmega \in \mathbb {R}^{p*N}\) is the analysis dictionary updated by GAP; \(\alpha \)/\(\alpha _m\) is sparse coefficients; \(\beta \)/\(\beta _m\) is the melioration for sparse coefficients based on nonlocal similarity and decided via the following scheme:

Step1: Find out Q similar image patches that have smallest Euclidean distance with the input patch in the whole image;

Step2: Calculate the similarity via formula (9) [8], define the set of similar blocks \(\theta '\):

$$\begin{aligned} \theta ':=\{j|w(u_{in},u_j^{le})>\varsigma \}, \end{aligned}$$
(8)
$$\begin{aligned} w(u_{in},u_j^{le})=\dfrac{1}{Z(m)}e^{-\dfrac{\left\| u_m-u_n\right\| _{2,a}^2}{h^2}},\end{aligned}$$
(9)
$$\begin{aligned} Z(m)=\sum \limits _{n}e^{-\dfrac{\left\| u_m-u_n\right\| _{2,a}^2}{h^2}}, \end{aligned}$$
(10)

where \(u_{in}\) means the input image to get \(\beta \); \(u_j^{le}\) represents the \(j_{th}\) patch in Step 1; \(\varsigma \) is the threshold; a is the standard deviation of the Gaussian kernel; h is the scaling parameter; Z(m) is the normalized coefficient;

Step 3: Calculate the similar image patch:

$$\begin{aligned} \tilde{u}=\sum w(u_{in},u_j^{le})*u_j^{le},j\in \theta ', \end{aligned}$$
(11)

Step 4: Calculate \(\beta \):

$$\begin{aligned} \tilde{u}=D\beta . \end{aligned}$$
(12)
Table 1. The proposed method SRMM.
Full size table

To solve formula (7), we divide it into two sub-problems based on alternative optimization:

Sub-problem 1: with u fixed, solve \(\alpha \):

$$\begin{aligned} \alpha =\arg \min \limits _{\alpha }T_1(\alpha )=\min \limits _{\alpha }\left\| u-D\alpha \right\| _2^2+\mu \left\| \alpha -\beta \right\| _1, \end{aligned}$$
(13)

Here we use the Iterative Re-weighted Least Squares Minimization (IRLSM) [21] for simplicity.

Sub-problem 2: with \(\alpha \) fixed, solve u:

$$\begin{aligned} u=\arg \min \limits _{u}T_2(u)=\min \limits _{u}\left\| u-D\alpha \right\| _2^2+\tau (\left\| u-s^i\right\| _2^2+\dfrac{2}{\eta _i}\left\| y-u\right\| _2^2)+\lambda \left\| \varOmega u\right\| _0. \end{aligned}$$
(14)

We choose GAP to solve (14) as it’s the IRLSM method in cosparsity’s framework.

The complete method called SRMM can be found in Table 1.

4 Experiments

In this section, experiments are displayed to demonstrate the effectiveness and superiority of the proposed method (SRMM). Firstly, experiments are conducted on commonly used natural images in Fig. 1(a)–(d). We add different degrees of mixed noise (\(\sigma =5\) and \(\sigma =25\)) to test the availability. Examples are shown in Fig. 1(e)–(f). Moreover, two micro-focus X-ray noisy images are used as the real-data experimental subjects, see Fig. 1(g)–(h). Several methods are selected as comparisons: (1) mixed noise removal: PURE-LET (P-L) [1], GAT+BM3D (G+B) [2], GAT+BLS-GSM (G+B-G) [3], DAR [17]; (2) sparse-based: K-SVD [7] and AK-SVD [14]; (3) multiple-constraint-based: CSR [10].

Fig. 1.
figure 1

Original images and noisy images.

Full size image

Our proposed denoising approach has seven vital paramete6rs to be set: the patch size, the number of the atoms in learning dictionary, the number of the cluster centers, the iterations number \(T_{iter}\) and three regularization parameters \(\mu \), \(\tau \) and \(\lambda \). We set the patch size \(7\times 7\) and cluster centers number 8. Each cluster has 200 atoms. The regularization terms \(\mu \), \(\tau \) and \(\lambda \) are set to be 1, 0.5 and 0.4, respectively. The maximum total iterations number \(T_{iter}\) is set 10. These parameters are selected empirically for higher performance in our experiment, but suitable variations to these parameters are acceptable depending on different image size and noise strength.

Due to the space limitation, we only show partial enlarged results of two natural images in Figs. 2 and 3. For micro-focus X-ray images, we demonstrate the edge-detection results in Figs. 4 and 5 as the micro-focus X-ray images are low-contrast and hard to compare directly. Objective indexes can be seen in Tables 2, 3 and 4. Full-Reference PSNR and MSSIM are used for natural images. And no reference MSR for smooth regions, LS for detail regions [17] and BRISQUE [22] are applied in micro-focus X-ray images since it’s unable to get the noise-free micro-focus X-ray images.

Fig. 2.
figure 2

Partial enlarged reconstruction results of Barbara.

Full size image
Fig. 3.
figure 3

Partial enlarged reconstruction results of Lena.

Full size image
Fig. 4.
figure 4

Partial enlarged edge-detection results (Canny) of Capacitance.

Full size image
Fig. 5.
figure 5

Edge-detection results (local variance) of Bubbles.

Full size image

For natural images, PURE-LET will not only blur details, but also bring artifacts that debase reconstruction quality. GAT+BLS-GSM suffers from over-smoothing, lots of details are lost. DAR suffers from details loss with slight noise remain. GAT+BM3D outperforms PURE-LET, GAT+BLS-GSM and DAR in visual effects. However, some details, especially weak edges are over-smoothed. K-SVD can preserve most details but there is noise left. Results of AK-SVD are not ideal for the reason that the edge transitions are not smooth. As the noise level especially the Poisson noise strength grows stronger, the noise in the non-smooth area is hard to be removed sufficiently. With two constraints, CSR reaches better results than the above methods. Nevertheless, some details are still lost. As for the proposed method, it can remove most mixed Gaussian-Poisson noise and in the meantime, preserve more edges and details.

For the reconstruction results of micro-focus X-ray images in Fig. 4, PURE-LET fails to reconstruct ‘Capacitance’. GAT+BM3D, DAR, K-SVD, AK-SVD and CSR suffer from over-smoothness. The details of pins as indicated in Fig. 4 are less. The noise is not well eliminated in GAT+BLS-GSM. The proposed method SRMM, however, removes most noise and preserves more details in the pins. For ‘Bubbles’ in Fig. 5, the proposed method and GAT+BM3D can reconstruct most details including small and weak bubbles in the left side. But GAT+BM3D’s result is smoother than SRMM. Other methods fail to reconstruct the shapes of the weak and small bubbles as much as SRMM does.

Table 2. Objective indexes of natural images with \(\sigma =5\) in mixed noise.
Full size table

Tables 2, 3 and 4 show the objective indexes of all experiments. In Table 2, SRMM achieves higher PSNR and MSSIM in most results. Also, for micro-focus X-ray images, SRMM gets lowest BRISQUE indicating best overall quality. And relatively higher MSR (more noise removed) and higher LS (more structures kept) mean that SRMM can attain a better balance between noise-removal and edge-preservation.

Table 3. Objective indexes of natural images with \(\sigma =25\) in mixed noise.
Full size table
Table 4. Objective indexes of micro-focus X-ray images.
Full size table

5 Conclusions

In this paper, we propose a new image sparse reconstruction method for mixed Gaussian-Poisson noise with multiple constraints. In our model, three priors, namely synthesis sparsity, analysis sparsity and nonlocal similarity, are combined as the sparse constraints to obtain better results. Among them, the synthesis sparsity and analysis sparsity are complementary and reconstruct noisy images by means of re-covering the sparsity. And nonlocal similarity explores the structure information that helps to reconstruct details. The objective function is solved via IRLSM and GAP alternatively. Experimental results show that the proposed method outperforms the contrast methods both in visual effects and indexes, which will improve reconstruction quality.