Adaptive Image Denoising by Targeted Databases

I Introduction

Image denoising is a classical signal recovery problem where the goal is to restore a clean image from its observations. Although image denoising has been studied for decades, the problem remains a fundamental one as it is the test bed for a variety of image processing tasks.

For example, in non-local means (NLM) , $\Phi$ is a weighted average of the reference patches, whereas in BM3D , $\Phi$ is a transform-shrinkage operation.

I-B Internal vs External Denoising

For any patch-based denoising algorithm, the denoising performance is intimately related to the reference patches $\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{k}$. Typically, there are two sources of these patches: the noisy image itself and an external database of patches. The former is known as internal denoising , whereas the latter is known as external denoising .

Internal denoising is practically more popular than external denoising because it is computationally less expensive. Moreover, internal denoising does not require a training stage, hence making it free of training bias. Furthermore, Glasner et al. showed that patches tend to recur within an image, e.g., at a different location, orientation, or scale. Thus searching for patches in the noisy image is often a plausible approach. However, on the downside, internal denoising often fails for rare patches — patches that seldom recur in an image. This phenomenon is known as “rare patch effect”, and is widely regarded as a bottleneck of internal denoising . There are some works attempting to alleviate the rare patch problem. However, the extent to which these methods can achieve is still limited.

External denoising is an alternative solution to internal denoising. Levin et al. showed that in the limit, the theoretical minimum mean squared error of denoising is achievable using an infinitely large external database. Recently, Chan et al. developed a computationally efficient sampling scheme to reduce the complexity and demonstrated practical usage of large databases. However, in most of the recent works on external denoising, e.g., , the databases used are generic. These databases, although large in volume, do not necessarily contain useful information to denoise the noisy image of interest. For example, it is clear that a database of natural images is not helpful to denoise a noisy portrait image.

I-C Adaptive Image Denoising

In this paper, we propose an adaptive image denoising algorithm using a targeted external database instead of a generic database. Here, a targeted database refers to a database that contains images relevant to the noisy image only. As will be illustrated in later parts of this paper, targeted external databases could be obtained in many practical scenarios, such as text images (e.g., newspapers and documents), human faces (under certain conditions), and images captured by multiview camera systems. Other possible scenarios include images of license plates, medical CT and MRI images, and images of landmarks.

The concept of using targeted external databases has been proposed in various occasions, e.g., . However, none of these methods are tailored for image denoising problems. The objective of this paper is to bridge the gap by addressing the following question:

(Q): Suppose we are given a targeted external database, how should we design a denoising algorithm which can maximally utilize the database?

Here, we assume that the reference patches $\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{k}$ are given. We emphasize that this assumption is application specific — for the examples we mentioned earlier (e.g., text, multiview, face, etc), the assumption is typically true because these images have relatively less variety in content.

When the reference patches are given, question (Q) may look trivial at the first glance because we can extend existing internal denoising algorithms in a brute-force way to handle external databases. For example, one can modify existing algorithms, e.g., , so that the patches are searched from a database instead of the noisy image. Likewise, one can also treat an external database as a “video” and feed the data to multi-image denoising algorithms, e.g., . However, the problem of these approaches is that the brute force modifications are heuristic. There is no theoretical guarantee of performance. This suggests that a straight-forward modification of existing methods does not solve question (Q), as the database is not maximally utilized.

An alternative response to question (Q) is to train a statistical prior of the targeted database, e.g., . The merit of this approach is that the performance often has theoretical guarantee because the denoising problem can now be formulated as a maximum a posteriori (MAP) estimation. However, the drawback is that many of these methods require a large number of training samples which is not always available in practice.

I-D Contributions and Organization

In view of the above seemingly easy yet challenging question, we introduced a new denoising algorithm using targeted external databases in . Compared to existing methods, the method proposed in achieves better performance and only requires a small number of external images. In this paper, we extend by offering the following new contributions:

Generalization of Existing Methods. We propose a generalized framework which encapsulates a number of denoising algorithms. In particular, we show (in Section III-B) that the proposed group sparsity minimization generalizes both fixed basis and PCA methods. We also show (in Section IV-B) that the proposed local Bayesian MSE solution is a generalization of many spectral operations in existing methods.

Improvement Strategies. We propose two improvement strategies for the generalized denoising framework. In Section III-D, we present a patch selection optimization to improve the patch search process. In Section IV-D, we present a soft-thresholding and a hard-thresholding method to improve the spectral coefficients learned by the algorithm.

Detailed Proofs. Proofs of the results in this paper and are presented in the Appendix.

The rest of the paper is organized as follows. After outlining the design framework in Section II, we present the above contributions in Section III – IV. Experimental results are discussed in Section V, and concluding remarks are given in Section VI.

II Optimal Linear Denoising Filter

The foundation of our proposed method is the classical optimal linear denoising filter design problem . In this section, we give a brief review of the design framework and highlight its limitations.

subject to the constraint that $\boldsymbol{U}$ is an orthonormal matrix.

The joint optimization (3) can be solved by noting the following Lemma.

The proof of Lemma 1 is given in . With Lemma 1, the denoised patch as a consequence of (3) is as follows.

The denoised patch $\boldsymbol{\widehat{p}}$ using the optimal $\boldsymbol{U}$ and $\boldsymbol{\Lambda}$ of (3) is

where $\boldsymbol{U}$ is any orthonormal matrix with the first column $\boldsymbol{u}_{1}=\boldsymbol{p}/\|\boldsymbol{p}\|_{2}$.

Lemma 2 states that if hypothetically we are given the ground truth $\boldsymbol{p}$, the optimal denoising process is to first project the noisy observation $\boldsymbol{q}$ onto the subspace spanned by $\boldsymbol{p}$, then perform a Wiener shrinkage $\|\boldsymbol{p}\|^{2}/(\|\boldsymbol{p}\|^{2}+\sigma^{2})$, and finally re-project the shrinkage coefficients to obtain the denoised estimate. However, since in reality we never have access to the ground truth $\boldsymbol{p}$, this optimal result is not achievable.

II-B Problem Statement

Since the oracle optimal filter is not achievable in practice, the question becomes whether it is possible to find a surrogate solution that does not require the ground truth $\boldsymbol{p}$.

To answer this question, it is helpful to separate the joint optimization (3) by first fixing $\boldsymbol{U}$ and minimize the MSE with respect to $\boldsymbol{\Lambda}$. In this case, one can show that (4) achieves the minimum when

in which the minimum MSE estimator is given by

where $\{\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{d}\}$ are the columns of $\boldsymbol{U}$.

Inspecting (6), we identify two parts of the problem:

Determine $\boldsymbol{U}$. The choice of $\boldsymbol{U}$ plays a critical role in the denoising performance. In literature, $\boldsymbol{U}$ are typically chosen as the FFT or the DCT bases . In , the PCA bases of various data matrices are proposed. However, the optimality of these bases is not fully understood.

Determine $\boldsymbol{\Lambda}$. Even if $\boldsymbol{U}$ is fixed, the optimal $\boldsymbol{\Lambda}$ in (5) still depends on the unknown ground truth $\boldsymbol{p}$. In , $\boldsymbol{\Lambda}$ is determined by hard-thresholding a stack of DCT coefficients or applying an empirical Wiener filter constructed from a first-pass estimate. In , $\boldsymbol{\Lambda}$ is formed by the PCA coefficients of a set of relevant noisy patches. Again, it is unclear which of these is optimal.

Motivated by the problems about $\boldsymbol{U}$ and $\boldsymbol{\Lambda}$, in the following two sections we present our proposed method for each of these problems. We discuss its relationship to prior works, and present ways to further improve it.

III Determine 𝑼𝑼\boldsymbol{U}

In this section, we present our proposed method to determine the basis matrix $\boldsymbol{U}$ and show that it is a generalization of a number of existing denoising algorithms. We also discuss ways to improve $\boldsymbol{U}$.

Given a noisy patch $\boldsymbol{q}$ and a targeted database $\{\boldsymbol{p}_{j}\}_{j=1}^{n}$, our first task is to fetch the $k$ most “relevant” patches. The patch selection is performed by measuring the similarity between $\boldsymbol{q}$ and each of $\{\boldsymbol{p}_{j}\}_{j=1}^{n}$, defined as

We note that (7) is equivalent to the standard $k$ nearest neighbors ($k$NN) search.

III-B Group Sparsity

Without loss of generality, we assume that the $k$NN returned by the above procedure are the first $k$ patches of the data, i.e., $\{\boldsymbol{p}_{j}\}_{j=1}^{k}$. Our goal now is to construct $\boldsymbol{U}$ from $\{\boldsymbol{p}_{j}\}_{j=1}^{k}$.

where $\boldsymbol{P}\overset{\text{def}}{=}[\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{k}]$. The equality constraint in (9) ensures that $\boldsymbol{U}$ is orthonormal. Thus, the solution of (9) is an orthonormal matrix $\boldsymbol{U}$ which maximizes the group sparsity of the data $\boldsymbol{P}$.

Interestingly, and surprisingly, the solution of (9) is indeed identical to the classical principal component analysis (PCA). The following lemma summarizes the observation.

where $\boldsymbol{S}$ is the corresponding eigenvalue matrix.

In practice, it is possible to improve the fidelity of the data matrix $\boldsymbol{P}$ by introducing a diagonal weight matrix

for some user tunable parameter $h$ and a normalization constant $Z\overset{\text{def}}{=}\boldsymbol{1}^{T}\boldsymbol{W}\boldsymbol{1}$. Consequently, we can define

Hence (10) becomes $[\boldsymbol{U},\boldsymbol{S}]=\mbox{eig}(\boldsymbol{P}\boldsymbol{W}\boldsymbol{P}^{T})$.

III-C Relationship to Prior Works

The fact that (10) is the solution to a group sparsity minimization problem allows us to understand the performance of a number of existing denoising algorithms to some extent.

It is perhaps a misconception that the underlying principle of BM3D is to enforce sparsity of the 3-dimensional data volume (which we shall call it a 3-way tensor). However, what BM3D enforces is the group sparsity of the slices of the tensor, not the sparsity of the tensor.

To see this, we note that the 3-dimensional transforms in BM3D are separable (e.g., DCT2 + Haar in its default setting). If the patches $\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{k}$ are sufficiently similar, the DCT2 coefficients will be similar in both magnitude and location By DCT2 location we meant the frequency of the DCT2 components.. Therefore, by fixing the frequency location of a DCT2 coefficient and tracing the DCT2 coefficients along the third axis, the output signal will be almost flat. Hence, the final Haar transform will return a sparse vector. Clearly, such sparsity is based on the stationarity of the DCT2 coefficients along the third axis. In essence, this is group sparsity.

III-C2 HOSVD [9]

where $\times_{k}$ denotes a tensor mode-$k$ multiplication .

As reported in , the performance of HOSVD is indeed worse than BM3D. This phenomenon can now be explained, because HOSVD ignores the fact that image patches tend to be group sparse instead of being tensor sparse.

III-C3 Shape-adaptive BM3D [4]

Consequently, the PCA of $\overline{\boldsymbol{P}}$ is equivalent to SA-BM3D. Here, the matrix $\boldsymbol{W}_{s}$ is used to control the relative emphasis of each pixel in the spatial coordinate.

III-C4 BM3D-PCA [5] and LPG-PCA [7]

The idea of both BM3D-PCA and LPG-PCA is that given $\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{k}$, $\boldsymbol{U}$ is determined as the principal components of $\boldsymbol{P}=[\boldsymbol{p}_{1},\ldots,\boldsymbol{p}_{k}]$. Incidentally, such approaches arrive at the same result as (10), i.e., the principal components are indeed the solution of a group sparse minimization. However, the key of using the group sparsity is not noticed in and . This provides additional theoretical justifications for both methods.

III-C5 KSVD [18]

In KSVD, the dictionary plays the role of our basis matrix $\boldsymbol{U}$. The dictionary can be trained either from the single noisy image, or from an external (generic or targeted) database. However, the training is performed once for all patches of the image. In other words, the noisy patches share a common dictionary. In our proposed method, each noisy patch has an individually trained basis matrix. Clearly, the latter approach, while computationally more expensive, is significantly more data adaptive than KSVD.

III-D Improvement: Patch Selection Refinement

The optimization problem (9) suggests that the $\boldsymbol{U}$ computed from (10) is the optimal basis with respect to the reference patches $\{\boldsymbol{p}_{j}\}_{j=1}^{k}$. However, one issue that remains is how to improve the selection of $k$ patches from the original $n$ patches. Our proposed approach is to formulate the patch selection as an optimization problem

where $\boldsymbol{c}=[c_{1},\cdots,c_{n}]^{T}$ with $c_{j}\overset{\text{def}}{=}\|\boldsymbol{q}-\boldsymbol{p}_{j}\|_{2}$, $\varphi(\boldsymbol{x})$ is a penalty function and $\tau>0$ is a parameter. In (14), each $c_{j}$ is the distance $\|\boldsymbol{q}-\boldsymbol{p}_{j}\|_{2}$, and $x_{j}$ is a weight indicating the emphasis of $\|\boldsymbol{q}-\boldsymbol{p}_{j}\|_{2}$. Therefore, the minimizer of (14) is a sequence of weights that minimize the overall distance.

To gain more insight into (14), we first consider the special case when the penalty term $\varphi(\boldsymbol{x})=0$. We claim that, under this special condition, the solution of (14) is equivalent to the original $k$NN solution in (7). This result is important, because $k$NN is a fundamental building block of all patch-based denoising algorithms. By linking $k$NN to the optimization formulation in (14) we provide a systematic strategy to improve the $k$NN.

The proof of the equivalence between $k$NN and (14) can be understood via the following case study where $n=2$ and $k=1$. In this case, the constraints $\boldsymbol{x}^{T}\boldsymbol{1}=1$ and $0\leq\boldsymbol{x}\leq 1$ form a closed line segment in the positive quadrant. Since the objective function $\boldsymbol{c}^{T}\boldsymbol{x}$ is linear, the optimal point must be at one of the vertices of the line segment, which is either $\boldsymbol{x}=^{T}$, or $\boldsymbol{x}=^{T}$. Thus, by checking which of $c_{1}$ or $c_{2}$ is smaller, we can determine the optimal solution by setting $x_{1}=1$ if $c_{1}$ is smaller (and vice versa). Correspondingly, if $x_{1}=1$, then the first patch $\boldsymbol{p}_{1}$ should be selected. Clearly, the solution returned by the optimization is exactly the $k$NN solution. A similar argument holds for higher dimensions, hence justifies our claim.

Knowing that $k$NN can be formulated as (14), our next task is to choose an appropriate penalty term. The following are two possible choices.

The penalized problem (15) suggests that the optimal $k$ reference patches should not be determined merely from $\|\boldsymbol{q}-\boldsymbol{p}_{j}\|_{2}$ (which could be problematic due to the noise present in $\boldsymbol{q}$). Instead, a good reference patch should also be similar to all other patches that are selected. The cross similarity term $x_{i}x_{j}\|\boldsymbol{p}_{i}-\boldsymbol{p}_{j}\|_{2}$ provides a way for such measure. This shares some similarities to the patch ordering concept proposed by Cohen and Elad . The difference is that the patch ordering proposed in is a shortest path problem that tries to organize the noisy patches, whereas ours is to solve a regularized optimization.

Problem (15) is in general not convex because the matrix $\boldsymbol{B}$ is not positive semidefinite. One way to relax the formulation is to consider $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$. Geometrically, the solution of using $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$ tends to identify patches that are close to the sum of all other patches in the set. In many cases, this is similar to $\varphi(\boldsymbol{x})=\boldsymbol{x}^{T}\boldsymbol{B}\boldsymbol{x}$ which finds patches that are similar to every individual patch in the set. In practice, we find that the difference between $\varphi(\boldsymbol{x})=\boldsymbol{x}^{T}\boldsymbol{B}\boldsymbol{x}$ and $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$ in the final denoising result (PSNR of the entire image) is marginal. Thus, for computational efficiency we choose $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$.

III-D2 Regularization by First-pass Estimate

The second choice of $\varphi(\boldsymbol{x})$ is based on a first-pass estimate $\overline{\boldsymbol{p}}$ using some denoising algorithms, for example, BM3D or the proposed method without this patch selection step. In this case, by defining $e_{j}\overset{\text{def}}{=}\|\overline{\boldsymbol{p}}-\boldsymbol{p}_{j}\|_{2}$ we consider the penalty function $\varphi(\boldsymbol{x})=\boldsymbol{e}^{T}\boldsymbol{x}$, where $\boldsymbol{e}=[e_{1},\cdots,e_{n}]^{T}$. This implies the following optimization problem

By identifying the objective of (16) as $(\boldsymbol{c}+\tau\boldsymbol{e})^{T}\boldsymbol{x}$, we observe that (16) can be solved in closed form by locating the $k$ smallest entries of the vector $\boldsymbol{c}+\tau\boldsymbol{e}$.

The interpretation of (16) is straight-forward: The linear combination of $\|\boldsymbol{q}-\boldsymbol{p}_{j}\|_{2}$ and $\|\overline{\boldsymbol{p}}-\boldsymbol{p}_{j}\|_{2}$ shows a competition between the noisy patch $\boldsymbol{q}$ and the first-pass estimate $\overline{\boldsymbol{p}}$. In most of the common scenarios, $\|\boldsymbol{q}-\boldsymbol{p}_{j}\|_{2}$ is preferred when noise level is low, whereas $\overline{\boldsymbol{p}}$ is preferred when noise level is high. This in turn requires a good choice of $\tau$. Empirically, we find that $\tau=0.01$ when $\sigma<30$ and $\tau=1$ when $\sigma>30$ is a good balance between the performance and generality.

III-D3 Comparisons

To demonstrate the effectiveness of the two proposed patch selection steps, we consider a ground truth (clean) patch shown in Figure 2 (a). From a pool of $n=200$ reference patches, we apply an exhaustive search algorithm to choose $k=40$ patches that best match with the noisy observation $\boldsymbol{q}$, where the first 10 patches are shown in Figure 2 (b). The results of the two selection refinement methods are shown in Figure 2 (c)-(d), where in both cases the parameter $\tau$ is adjusted for the best performance. For the case of $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$, we set $\tau=1/(200n)$ when $\sigma<30$ and $\tau=1/(2n)$ when $\sigma>30$. For the case of $\varphi(\boldsymbol{x})=\boldsymbol{e}^{T}\boldsymbol{x}$, we use the denoised result of BM3D as the first-pass estimate $\overline{\boldsymbol{p}}$, and set $\tau=0.01$ when $\sigma<30$ and $\tau=1$ when $\sigma>30$. The results in Figure 3 show that the PSNR increases from 28.29 dB to 28.50 dB if we use $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$, and further increases to 29.30 dB if we use $\varphi(\boldsymbol{x})=\boldsymbol{e}^{T}\boldsymbol{x}$. The full performance comparison is shown in Figure 4, where we show the PSNR curve for a range of noise levels of an image. Since the performance of $\varphi(\boldsymbol{x})=\boldsymbol{e}^{T}\boldsymbol{x}$ is consistently better than $\varphi(\boldsymbol{x})=\boldsymbol{1}^{T}\boldsymbol{B}\boldsymbol{x}$, in the rest of the paper we focus on $\varphi(\boldsymbol{x})=\boldsymbol{e}^{T}\boldsymbol{x}$.

IV Determine 𝚲𝚲\boldsymbol{\Lambda}

In this section, we present our proposed method to determine $\boldsymbol{\Lambda}$ for a fixed $\boldsymbol{U}$. Our proposed method is based on the concept of a Bayesian MSE estimator.

Assuming that the prior distribution $f(\boldsymbol{p})$ is known, it is natural to consider the Bayesian mean squared error (BMSE) between the estimate $\boldsymbol{\widehat{p}}\overset{\text{def}}{=}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{T}\boldsymbol{q}$ and the ground truth $\boldsymbol{p}$:

The BMSE defined in (18) suggests that the optimal $\boldsymbol{\Lambda}$ should be the minimizer of the optimization problem

In the next subsection we discuss how to solve (19).

IV-B Localized Prior from the Targeted Database

Minimizing BMSE over $\boldsymbol{\Lambda}$ involves knowing the prior distribution $f(\boldsymbol{p})$. However, in general, the exact form of $f(\boldsymbol{p})$ is never known. This leads to many popular models in the literature, e.g., Gaussian mixture model , the field of expert model , and the expected patch log-likelihood model (EPLL) .

One common issue of all these models is that the prior $f(\boldsymbol{p})$ is built from a generic database of patches. In other words, the $f(\boldsymbol{p})$ models all patches in the database. As a result, $f(\boldsymbol{p})$ is often a high dimensional distribution with complicated shapes.

In our targeted database setting, the difficult prior modeling becomes a much simpler task. The reason is that while the shape of the distribution $f(\boldsymbol{p})$ is still unknown, the subsampled reference patches (which are few but highly representative) could be well approximated as samples drawn from a single Gaussian centered around some mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$. Therefore, by appropriately estimating $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ of this localized prior, we can derive the optimal $\boldsymbol{\Lambda}$ as given by the following Lemma:

Let $f(\boldsymbol{q}\,|\,\boldsymbol{p})=\mathcal{N}(\boldsymbol{p},\sigma^{2}\boldsymbol{I})$, and let $f(\boldsymbol{p})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ for any vector $\boldsymbol{\mu}$ and matrix $\boldsymbol{\Sigma}$, then the optimal $\boldsymbol{\Lambda}$ that minimizes (18) is

where $\boldsymbol{G}\overset{\text{def}}{=}\boldsymbol{U}^{T}\boldsymbol{\mu}\boldsymbol{\mu}^{T}\boldsymbol{U}+\boldsymbol{U}^{T}\boldsymbol{\Sigma}\boldsymbol{U}$.

To specify $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, we let

where $w_{j}$ is the $j$th diagonal entry of $\boldsymbol{W}$ defined in (11). Intuitively, an interpretation of (21) is that $\boldsymbol{\mu}$ is the non-local mean of the reference patches. However, the more important part of (21) is $\boldsymbol{\Sigma}$, which measures the uncertainty of the reference patches with respect to $\boldsymbol{\mu}$. This uncertainty measure makes some fundamental improvements to existing methods which will be discussed in Section IV-C.

We note that Lemma 4 holds even if $f(\boldsymbol{p})$ is not Gaussian. In fact, for any distribution $f(\boldsymbol{p})$ with the first cumulant $\boldsymbol{\mu}$ and the second cumulant $\boldsymbol{\Sigma}$, the optimal solution in (41) still holds. This result is equivalent to the classical linear minimum MSE (LMMSE) estimation .

From a computational perspective, $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ defined in (21) lead to a very efficient implementation as illustrated by the following lemma.

Using $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ defined in (21), the optimal $\boldsymbol{\Lambda}$ is given by

where $\boldsymbol{S}$ is the eigenvalue matrix of $\boldsymbol{P}\boldsymbol{W}\boldsymbol{P}^{T}$.

Combining Lemma 5 with Lemma 3, we observe that for any set of reference patches $\{\boldsymbol{p}_{j}\}_{j=1}^{k}$, $\boldsymbol{U}$ and $\boldsymbol{\Lambda}$ can be determined simultaneously through the eigen-decomposition of $\boldsymbol{P}\boldsymbol{W}\boldsymbol{P}^{T}$. Therefore, we arrive at the overall algorithm shown in Algorithm 1.

IV-C Relationship to Prior Works

BM3D and its variants have two denoising steps. In the first step, the algorithm applies a basis matrix $\boldsymbol{U}$ (either a pre-defined basis such as DCT, or a basis learned from PCA). Then, it applies a hard-thresholding to the projected coefficients to obtain a filtered image $\overline{\boldsymbol{p}}$. In the second step, the filtered image $\overline{\boldsymbol{p}}$ is used as a pilot estimate to the desired spectral component

Following our proposed Bayesian framework, we observe that the role of using $\overline{\boldsymbol{p}}$ in (23) is equivalent to assuming a dirac delta prior

In other words, the prior that BM3D assumes is concentrated at one point, $\overline{\boldsymbol{p}}$, and there is no measure of uncertainty. As a result, the algorithm becomes highly sensitive to the first-pass estimate. In contrast, (21) suggests that the first-pass estimate can be defined as a non-local mean solution. Additionally, we incorporate a covariance matrix $\boldsymbol{\Sigma}$ to measure the uncertainty of observing $\boldsymbol{\mu}$. These provide a more robust estimate to the denoising algorithm which is absent from BM3D and its variants.

IV-C2 LPG-PCA [7]

In LPG-PCA, the $i$th spectral component $\lambda_{i}$ is defined as

where $\boldsymbol{q}$ is the noisy patch. The (implicit) assumption in is that $(\boldsymbol{u}_{i}^{T}\boldsymbol{q})^{2}\approx(\boldsymbol{u}_{i}^{T}\boldsymbol{p})^{2}+\sigma^{2}$, and so by substituting $(\boldsymbol{u}_{i}^{T}\boldsymbol{p})^{2}\approx(\boldsymbol{u}_{i}^{T}\boldsymbol{q})^{2}-\sigma^{2}$ into (5) yields (25). However, the assumption implies the existence of a perturbation $\Delta\boldsymbol{p}$ such that $(\boldsymbol{u}_{i}^{T}\boldsymbol{q})^{2}=(\boldsymbol{u}_{i}^{T}(\boldsymbol{p}+\Delta\boldsymbol{p}))^{2}+\sigma^{2}$. Letting $\overline{\boldsymbol{p}}=\boldsymbol{p}+\Delta\boldsymbol{p}$, we see that LPG-PCA implicitly assumes a dirac prior as in (23) and (24). The denoising result depends on the magnitude of $\Delta\boldsymbol{p}$.

IV-C3 Generic Global Prior [22]

As a comparison to methods using generic databases such as , we note that the key difference lies in the usage of a global prior versus a local prior. Figure 5 illustrates the concept pictorially. The generic (global) prior $f(\boldsymbol{p})$ covers the entire space, whereas the targeted (local) prior is concentrated at its mean. The advantage of the local prior is that it allows one to denoise an image with few reference patches. It saves us from the intractable computation of learning the global prior, which is a high-dimensional non-parametric function.

IV-C4 Generic Local Prior – EPLL [19], K-SVD [35, 18]

Compared to learning-based methods that use local priors, such as EPLL and K-SVD , the most important merit of the proposed method is that it requires significantly fewer training samples. A thorough justification will be discussed in Section V.

IV-C5 PLOW [48]

PLOW has a similar design process as ours by considering the optimal filter. The major difference is that in PLOW, the denoising filter is derived from the full covariance matrices of the data and noise. As we will see in the next subsection, the linear denoising filter of our work is a truncated SVD matrix computed from a set of similar patches. The merit of the truncation is that it often reduces MSE in the bias-variance trade off .

IV-D Improving 𝚲𝚲\boldsymbol{\Lambda}

where $\gamma>0$ is the penalty parameter, and $\alpha\in\{0,\,1\}$ controls which norm to be used. The solution to the minimization of (26) is given by the following lemma.

V Experimental Results

In this section, we present a set of experimental results.

The methods we choose for comparison are BM3D , BM3D-PCA , LPG-PCA , NLM , EPLL and KSVD . Except for EPLL and KSVD, all other four methods are internal denoising methods. We re-implement and modify the internal methods so that patch search is performed over the targeted external databases. These methods are iterated for two times where the solution of the first step is used as a basic estimate for the second step. The specific settings of each algorithm are as follows:

BM3D : As a benchmark of internal denoising, we run the original BM3D code provided by the authorhttp://www.cs.tut.fi/~foi/GCF-BM3D/. Default parameters are used in the experiments, e.g., the search window is $39\times 39$. We have included a discussion in Section V-B about the influence of different search window size to the denoising performance. As for external denoising, we implement an external version of BM3D. To ensure a fair comparison, we set the search window identical to other external denoising methods.

BM3D-PCA and LPG-PCA : $\boldsymbol{U}$ is learned from the best $k$ external patches, which is the same as in our proposed method. $\boldsymbol{\Lambda}$ is computed following (23) for BM3D-PCA and (25) for LPG-PCA. In BM3D-PCA’s first step, the threshold is set to $2.7\sigma$.

EPLL : In EPLL, the default patch prior is learned from a generic database (200,000 $8\times 8$ patches). For a fair comparison, we train the prior distribution from our targeted databases using the same EM algorithm mentioned in .

KSVD : In KSVD, two dictionaries are trained including a global dictionary and a targeted dictionary. The global dictionary is trained from a generic database of 100,000 $8\times 8$ patches by the KSVD authors. The targeted dictionary is trained from a targeted database of 100,000 $8\times 8$ patches containing similar content of the noisy image. Both dictionaries are of size $64\times 256$.

To emphasize the difference between the original algorithms (which are single-image denoising algorithms) and the corresponding new implementations for external databases, we denote the original, (single-image) denoising algorithms with “$i$” (internal), and the corresponding new implementations for external databases with “$e$” (external).

We add zero-mean Gaussian noise with standard deviations from $\sigma=20$ to $\sigma=80$ to the test images. The patch size is set as $8\times 8$ ($i.e.,d=64$), and the sliding step size is 6 in the first step and 4 in the second step. Two quality metrics, namely Peak Signal to Noise Ratio (PSNR) and Structural Similarity (SSIM) are used to evaluate the objective quality of the denoised images.

V-B Denoising Text and Documents

Our first experiment considers denoising a text image. The purpose is to simulate the case where we want to denoise a noisy document with the help of other similar but non-identical texts. This idea can be easily generalized to other scenarios such as handwritten signatures, bar codes and license plates.

To prepare this scenario, we capture randomly 8 regions of a document and add noise. We then build the targeted external database by cropping 9 arbitrary portions from a different document but with the same font sizes.

Figure 7 shows the denoising results when we add excessive noise ($\sigma=100$) to one query image. Among all the methods, the proposed method yields the highest PSNR and SSIM values. The PSNR is 5 dB better than the benchmark BM3D (internal) denoising algorithm. Some existing learning-based methods, such as EPLL, do not perform well due to the insufficient training samples from the targeted database. Compared to other external denoising methods, the proposed method shows a better utilization of the targeted database.

Since the default search window size for internal BM3D is only $39\times 39$, we further conduct experiments to explore the effect of different search window sizes for BM3D. The PSNR results are shown in Table I. We see that a larger window size improves the BM3D denoising performance since more patch redundancy can be exploited. However, even if we extend the search to an external database (which is the case for eBM3D), the performance is still worse than the proposed method.

In Figure 8, we plot and compare the average PSNR values on 8 test images over a range of noise levels. We observe that at low noise levels ($\sigma<30$), our proposed method performs worse than eBM3D-PCA and eLPG-PCA. One reason is that the patch variety of the text image database makes our estimate of $\boldsymbol{\Lambda}$ in $\eqref{eq:solution for mLambda}$ worse than the other two estimates in (23) and (25). However, as noise level increases, our proposed method outperforms other methods, which suggests that the prior of our method is more informative. For example, for $\sigma=60$, our average PSNR result is 1.26 dB better than the second best result by eBM3D-PCA.

For the two learning-based methods, i.e., EPLL and KSVD, as can be seen, using a targeted database yields better results than using a generic database, which validates the usefulness of a targeted database. However, they perform worse than other non-learning methods. One reason is that a large number of training samples are needed for these learning-based methods – for EPLL, the large number of samples is needed to build the Gaussian mixtures, whereas for KSVD, the large number of samples is needed to train the over-complete dictionary. In contrast, the proposed method is fully functional even if the database is small.

V-B2 Database Quality

The average patch-database distance is then defined as $\overline{d}(\mathcal{P})\overset{\text{def}}{=}(1/m)\sum_{i=1}^{m}d(\boldsymbol{p}_{i},\mathcal{P})$. Therefore, a smaller $\overline{d}(\mathcal{P})$ indicates that the database is more relevant to the ground truth (clean) image.

Figure 9 shows the results of six databases $\mathcal{P}$, where each is a random subset of the original targeted database. For all noise levels ($\sigma=$ 20 to 80), PSNR decreases linearly as the patch-to-database distance increase, Moreover, the decay rate is slower for higher noise levels. The result suggests that the quality of the database has a more significant impact under low noise conditions, and less under high noise conditions.

V-C Denoising Multiview Images

Our second experiment considers the scenario of capturing images using a multiview camera system. The multiview images are captured at different viewing positions. Suppose that one or more cameras are not functioning properly so that some images are corrupted with noise. Our goal is to demonstrate that with the help of the other clean views, the noisy view could be restored.

To simulate the experiment, we download 4 multivew datasets from Middlebury Computer Vision Pagehttp://vision.middlebury.edu/stereo/. Each set of images consists of 5 views. We add i.i.d. Gaussian noise to one view and then use the rest 4 views to assist in denoising.

In Figure 10, we visually show the denoising results of the “Barn” and “Cone” multiview datasets. In comparison to the competing methods, our proposed method has the highest PSNR values. The magnified areas indicate that our proposed method removes the noise significantly and better reconstructs some fine details. In Figure 11, we plot and compare the average PSNR values on 4 test images over a range of noise levels. The proposed method is consistently better than its competitors. For example, for $\sigma=50$, our proposed method is 1.06 dB better than eBM3D-PCA and 2.73 dB better than iBM3D. The superior performance confirms our belief that with a good database, not any denoising algorithm would perform equally well. In fact, we still have to carefully design the denoising algorithm in order to maximize the performance by fully utilizing the database.

V-D Denoising Human Faces

Our third experiment considers denoising human face images. In low light conditions, images captured are typically corrupted by noise. To facilitate other high-level vision tasks such as recognition and tracking, denoising is a necessary pre-processing step. This experiment demonstrates the ability of denoising face images.

In this experiment, we use the Gore face database from , of which some examples are shown in the top row of Figure 12 (each image is $60\times 80$). We simulate the denoising task by adding noise to 8 randomly chosen images and then use the other images (29 other face images in our experiment) in the database to assist in denoising.

In the top row of Figure 12, we show some clean face images in the database while in the bottom row, we show one of the noisy faces and its denoising results (magnified). We observe that while the facial expressions are different and there are misalignments between images, the proposed method still generates robust results. In Figure 13, we plot the average PSNR curves on the 8 test images, where we see consistent gain compared to other methods.

V-E Runtime Comparison

Our current implementation is in MATLAB (single thread). The runtime is about 144s to denoise an image ($301\times 218$) with a targeted database consisting of 9 images of similar sizes. The code is run on an Intel Core i7-3770 CPU. In Table II, we show a runtime comparison with other methods. We observe that the runtime of the proposed method is indeed not significantly worse than other external methods. In particular, the runtime of the proposed method is in the same order of magnitude as eNLM, eBM3D, eBM3D-PCA and eLPG-PCA.

We note that most of the runtime of the proposed method is spent on searching similar patches and computing SVD. Speed improvement for the proposed method is possible. First, we can apply techniques to enable fast patch search, e.g., patch match , KD tree , or fast SVD . Second, random sampling schemes can be applied to further reduce the computational complexity . Third, since the denoising is independently performed on each patch, GPU can be used to parallelize the computation.

V-F Discussion and Future Work

One important aspect of the algorithm that we did not discuss in depth is the sensitivity. In particular, two questions must be answered. First, assuming that there is a perturbation on the database, how much MSE will be changed? Answering the question will provide us information about the sensitivity of the algorithm when there are changes in font size (in the text example), view angle (in the multiview example), and facial expression (in the face example). Second, given a clean patch, how many patches do we need to put in the database in order to ensure that the clean patch is close to at least one of the patches in the database? The answer to this question will inform us about the size of the targeted database. Both problems will be studied in our future work.

VI Conclusion

Classical image denoising methods based on a single noisy input or generic databases are approaching their performance limits. We proposed an adaptive image denoising algorithm using targeted databases. The proposed method applies a group sparsity minimization and a localized prior to learn the basis matrix and the spectral coefficients of the optimal denoising filter, respectively. We show that the new method generalizes a number of existing patch-based denoising algorithms such as BM3D, BM3D-PCA, Shape-adaptive BM3D, LPG-PCA, and EPLL. Based on the new framework, we proposed improvement schemes, namely an improved patch selection procedure for determining the basis matrix and a penalized minimization for determining the spectral coefficients. For a variety of scenarios including text, multiview images and faces, we demonstrated empirically that the proposed method has superior performance over existing methods. With the increasing amount of image data available online, we anticipate that the proposed method is an important first step towards a data-dependent generation of denoising algorithms.

Appendix A Appendix

Since each term in the sum of the objective function is non-negative, we can consider the minimization over each individual term separately. This gives

In (29), we temporarily dropped the orthogonality constraint $\boldsymbol{u}_{i}^{T}\boldsymbol{u}_{j}=0$, which will be taken into account later. The Lagrangian function of (29) is

where $\beta$ is the Lagrange multiplier. Differentiating $\mathcal{L}$ with respect to $\boldsymbol{u}_{i}$, $\lambda_{i}$ and $\beta$ yields

Setting $\partial\mathcal{L}/\partial\lambda_{i}=0$ yields

Substituting this $\lambda_{i}$ into (31) and setting $\partial\mathcal{L}/\partial\boldsymbol{u}_{i}=0$ yields

Therefore, the optimal pair ($\boldsymbol{u}_{i}$, $\beta$) of (29) must be the solution of (34). The corresponding $\lambda_{i}$ can be calculated via (33).

Referring to (34), we observe two possible scenarios. First, if $\boldsymbol{u}_{i}$ is any unit vector orthogonal to $\boldsymbol{p}_{i}$, and $\beta=0$, then (34) can be satisfied. This is a trivial solution, because $\boldsymbol{u}_{i}\bot\boldsymbol{p}$ implies $\boldsymbol{u}_{i}^{T}\boldsymbol{p}=0$, and hence $\lambda_{i}=0$. The second case is that

Substituting (35) shows that (34) is satisfied. This is the non-trivial solution. The corresponding $\lambda_{i}$ in this case is $\|\boldsymbol{p}\|^{2}/(\|\boldsymbol{p}\|^{2}+\sigma^{2})$.

Finally, taking into account of the orthogonality constraint $\boldsymbol{u}_{i}^{T}\boldsymbol{u}_{j}=0$ if $i\not=j$, we can choose $\boldsymbol{u}_{1}=\boldsymbol{p}/\|\boldsymbol{p}\|_{2}$, and $\boldsymbol{u}_{2}\bot\boldsymbol{u}_{1}$, $\boldsymbol{u}_{3}\bot\{\boldsymbol{u}_{1},\boldsymbol{u}_{2}\}$, $\ldots$, $\boldsymbol{u}_{d}\bot\{\boldsymbol{u}_{1},\boldsymbol{u}_{2},\ldots\boldsymbol{u}_{d-1}\}$. Therefore, the denoising result is

where $\boldsymbol{U}$ is any orthonormal matrix with the first column $\boldsymbol{u}_{1}=\boldsymbol{p}/\|\boldsymbol{p}\|_{2}$. ∎

A-B Proof of Lemma 3

Let $\boldsymbol{u}_{i}$ be the $i$th column of $\boldsymbol{U}$. Then, (9) becomes

Since each term in the sum of (36) is non-negative, we can consider each individual term

The constrained problem (37) can be solved by considering the Lagrange function,

Taking derivatives $\frac{\partial\mathcal{L}}{\partial\boldsymbol{u}_{i}}=0$ and $\frac{\partial\mathcal{L}}{\partial\beta}=0$ yield

Therefore, $\boldsymbol{u}_{i}$ is the eigenvector of $\boldsymbol{P}\boldsymbol{P}^{T}$, and $\beta$ is the corresponding eigenvalue. Since the eigenvectors are orthonormal to each other, the solution automatically satisfies the orthogonality constraint that $\boldsymbol{u}_{i}^{T}\boldsymbol{u}_{j}=0$ if $i\not=j$. ∎

A-C Proof of Lemma 4

First, by plugging $\boldsymbol{q}=\boldsymbol{p}+\boldsymbol{\eta}$ into BMSE we get

where $\boldsymbol{G}\overset{\text{def}}{=}\boldsymbol{U}^{T}\boldsymbol{\mu}\boldsymbol{\mu}^{T}\boldsymbol{U}+\boldsymbol{U}^{T}\boldsymbol{\Sigma}\boldsymbol{U}$ and $g_{i}$ is the $i$th diagonal entry in $\boldsymbol{G}$.

Therefore, the optimal $\lambda_{i}$ is $g_{i}/(g_{i}+\sigma^{2})$ and the optimal $\boldsymbol{\Lambda}$ is

A-D Proof of Lemma 5

First, we write $\boldsymbol{\Sigma}$ in (21) in the matrix form

It is not difficult to see that $\boldsymbol{1}^{T}\boldsymbol{W}\boldsymbol{P}^{T}=\boldsymbol{\mu}^{T},\boldsymbol{P}\boldsymbol{W}\boldsymbol{1}=\boldsymbol{\mu}$ and $\boldsymbol{1}^{T}\boldsymbol{W}\boldsymbol{1}=1$. Therefore,

Note that $\boldsymbol{G}=\boldsymbol{U}^{T}\boldsymbol{\mu}\boldsymbol{\mu}^{T}\boldsymbol{U}+\boldsymbol{U}^{T}\boldsymbol{\Sigma}\boldsymbol{U}=\boldsymbol{U}^{T}(\boldsymbol{\mu}\boldsymbol{\mu}^{T}+\boldsymbol{\Sigma})\boldsymbol{U}$. Substituting (42) into $\boldsymbol{G}$ and using equation (10), we have

A-E Proof of Lemma 6

Therefore, the minimization of (26) becomes

where $\gamma\|\Lambda\boldsymbol{1}\|_{\alpha}=\gamma\sum_{i=1}^{d}|\lambda_{i}|$ or $\gamma\sum_{i=1}^{d}\mathds{1}(\lambda_{i}\neq 0)$ for $\alpha=1$ or . We note that when $\alpha=$ 1 or 0, (44) is the standard shrinkage problem , in which a closed form solution exists. Following from , the solutions are given by

References