External Prior Guided Internal Prior Learning for Real-World Noisy Image Denoising

TL;DR

This paper introduces a novel denoising approach that combines external clean image priors with internal image priors learned from the noisy image itself, effectively handling complex real-world noise.

Contribution

It proposes an external prior guided internal prior learning framework that adaptively combines external and internal priors for improved real-world noisy image denoising.

Findings

Outperforms state-of-the-art denoising methods on real-world datasets.

Effectively models complex real-world noise beyond simple distributions.

Demonstrates the benefit of combining external and internal priors.

Abstract

Most of existing image denoising methods learn image priors from either external data or the noisy image itself to remove noise. However, priors learned from external data may not be adaptive to the image to be denoised, while priors learned from the given noisy image may not be accurate due to the interference of corrupted noise. Meanwhile, the noise in real-world noisy images is very complex, which is hard to be described by simple distributions such as Gaussian distribution, making real-world noisy image denoising a very challenging problem. We propose to exploit the information in both external data and the given noisy image, and develop an external prior guided internal prior learning method for real-world noisy image denoising. We first learn external priors from an independent set of clean natural images. With the aid of learned external priors, we then learn internal priors from the given noisy image to refine the prior model. The external and internal priors are formulated as a set of orthogonal dictionaries to efficiently reconstruct the desired image. Extensive experiments are performed on several real-world noisy image datasets. The proposed method demonstrates highly competitive denoising performance, outperforming state-of-the-art denoising methods including those designed for real-world noisy images.

Figures

Tables

Algorithm 1: External Prior Guided Internal Prior Learning

Input: Matrices ${\overline{𝒀}}_n$, external sub-dictionary ${𝑫}_{\text{E}}$, parameter vector $𝝀$

Initialization: initialize ${𝑫}^{\left(0\right)}={𝑼}_k$ by Eq. (2);

for $t=0,1,\dots ,T-1$ do

1. Update ${𝜶}_{n,m}^{\left(t+1\right)}$ by Eq. (7);

2. Update ${𝑫}_{\text{I}}^{\left(t+1\right)}$ by Eq. (9);

end for

Output: Internal orthogonal dictionary ${𝑫}_{\text{I}}^{\left(T\right)}$ and sparse codes ${\text{A}}^{\left(T\right)}$.

Algorithm 2: External Prior Guided Internal Prior Learning for

Real-World Noisy Image Denoising

Input: Noisy image $𝐲$, external PG prior GMM model

Initialization: ${\widehat{𝐱}}^{\left(0\right)}=𝐲$;

for $Ite=1:IteNum$ do

1. Extracting internal PGs ${\left\{{𝒀}_n\right\}}_{n=1}^N$ from ${\widehat{𝐱}}^{\left(Ite-1\right)}$;

Guided Internal Subspace Clustering:

for each PG ${𝒀}_n$ do

2. Calculate group mean ${𝝁}_n$ and form mean subtracted PG ${\overline{𝒀}}_n$;

3. Subspace clustering via Eq. (3);

end for

Guided Internal Orthogonal Dictionary Learning:

for the PGs in each subspace do

4. External PG prior guided internal orthogonal dictionary learning by

solving Eq. (5);

5. Recover each patch in all PGs via Eq. (13);

end for

6. Aggregate the recovered PGs of all subspaces to form the recovered

image ${\widehat{𝐱}}^{\left(Ite\right)}$;

end for

Output: The denoised image $\widehat{𝐱}$.

Table 3. Table I: Average PSNR (dB) and Run Time (seconds) of the “External”, “Internal”, and “Guided Internal” methods on 60 real-world noisy images (of size 500 × 500 × 3 500 500 3 500\times 500\times 3 ) cropped from Dataset 1 [ 42 ] .

	Noisy	External	Internal	Guided Internal
PSNR	34.51	38.21	38.07	38.75
Time	—	21.19	312.67	22.26

Table 4. Table II: PSNR(dB) results and Speed (sec.) of different methods on 15 cropped real-world noisy images used in [ 42 ] .

Setting	GAT-BM3D	CBM3D	WNNM	TID	MLP	CSF	TNRD	DnCNN	NI	NC	CC	Ours
Canon 5D	31.23	39.76	37.51	37.22	39.00	35.68	39.51	37.26	37.68	38.76	38.37	40.50
\cdashline2-13 ISO = 3200	30.55	36.40	33.86	34.54	36.34	34.03	36.47	34.13	34.87	35.69	35.37	37.05
\cdashline2-13	27.74	36.37	31.43	34.25	36.33	32.63	36.45	34.09	34.77	35.54	34.91	36.11
Nikon D600	28.55	34.18	33.46	32.99	34.70	31.78	34.79	33.62	34.12	35.57	34.98	34.88
\cdashline2-13 ISO = 3200	32.01	35.07	36.09	34.20	36.20	35.16	36.37	34.48	35.36	36.70	35.95	36.31
\cdashline2-13	39.78	37.13	39.86	35.58	39.33	39.98	39.49	35.41	38.68	39.28	41.15	39.23
Nikon D800	32.24	36.81	36.35	34.94	37.95	34.84	38.11	35.79	37.34	38.01	37.99	38.40
\cdashline2-13 ISO = 1600	33.86	37.76	39.99	35.19	40.23	38.42	40.52	36.08	38.57	39.05	40.36	40.92
\cdashline2-13	33.90	37.51	37.15	35.26	37.94	35.79	38.17	35.48	37.87	38.20	38.30	38.97
Nikon D800	36.49	35.05	38.60	33.70	37.55	38.36	37.69	34.08	36.95	38.07	39.01	38.66
\cdashline2-13 ISO = 3200	32.91	34.07	36.04	31.04	35.91	35.53	35.90	33.70	35.09	35.72	36.75	37.07
\cdashline2-13	40.20	34.42	39.73	33.07	38.15	40.05	38.21	33.31	36.91	36.76	39.06	38.52
Nikon D800	29.84	31.13	33.29	29.40	32.69	34.08	32.81	29.83	31.28	33.49	34.61	33.76
\cdashline2-13 ISO = 6400	27.94	31.22	31.16	29.86	32.33	32.13	32.33	30.55	31.38	32.79	33.21	33.43
\cdashline2-13	29.15	30.97	31.98	29.21	32.29	31.52	32.29	30.09	31.40	32.86	33.22	33.58
Average	32.43	35.19	35.77	33.36	36.46	35.33	36.61	33.86	35.49	36.43	36.88	37.15
Time (s)	10.9	6.9	151.5	7353.2	16.8	19.3	5.1	79.2	0.6	15.3	NA	23.9

Table 5. Table III: Average PSNR(dB) results of different methods on 60 real-world noisy images cropped from [ 42 ] .

Methods	GAT-BM3D	CBM3D	WNNM	MLP	CSF	TNRD	DnCNN	NI	NC	Ours
PSNR	34.33	36.34	37.67	38.13	37.40	38.32	34.99	36.53	37.57	38.75

Table 6. Table IV: Average PSNR(dB) results of different methods on the 1000 real-world noisy images from the DND dataset [ 56 ] .

Methods	GAT-BM3D	CBM3D	WNNM	MLP	CSF	TNRD	DnCNN	NI	NC	Ours
PSNR	30.07	32.14	33.28	34.02	33.87	34.15	32.41	35.11	36.07	36.41

Table 7. Table V: Average PSNR(dB) results of different methods on 100 real-world noisy images cropped from our new dataset.

Methods	GAT-BM3D	CBM3D	WNNM	MLP	CSF	TNRD	DnCNN	NI	NC	Ours
PSNR	33.54	37.14	35.18	37.34	37.07	37.48	34.74	35.70	36.76	37.64

Table 8. Table VI: Average Speed (sec.) results of different methods on 100 real-world noisy images cropped from our new dataset.

Methods	GAT-BM3D	CBM3D	WNNM	MLP	CSF	TNRD	DnCNN	NI	NC	Ours
Time	11.1	6.9	152.2	17.1	19.5	5.2	79.5	0.6	15.6	24.1

Equations

$$\text{vspace}-2mm\begin{array}{r}\ln \mathcal{L}=\sum _{l=1}^L\ln (\sum _{k=1}^K{\pi }_k\prod _{m=1}^M\mathcal{N}({\overline{\mathbf{x}}}_{l,m}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)).\end{array}$$

$$\text{vspace}-2mm{\mathbf{\Sigma }}_k={\mathbf{U}}_k{\mathbf{S}}_k{\mathbf{U}}_k^{\top }.$$

$$\text{vspace}-2mmP(k|{\overline{\mathbf{Y}}}_n)=\frac{{\prod }_{m=1}^M\mathcal{N}({\overline{\mathbf{y}}}_{n,m}|0,{\mathbf{\Sigma }}_k)}{{\sum }_{l=1}^K{\prod }_{m=1}^M\mathcal{N}({\overline{\mathbf{y}}}_{n,m}|0,{\mathbf{\Sigma }}_l)}$$

$$\text{vspace}-2mm{\mathbf{D}}_k\triangleq [{\mathbf{D}}_{k,\text{E}}{\mathbf{D}}_{k,\text{I}}]\in {\mathbb{R}}^{3p^2\times 3p^2},$$

$$\text{vspace}-2mm\begin{array}{rl}\underset{{\mathbf{D}}_{\text{I}},\{{\mathbf{\alpha }}_{n,m}\}}{\min }& \sum _{n=1}^N\sum _{m=1}^M(\Vert {\overline{\mathbf{y}}}_{n,m}-\mathbf{D}{\mathbf{\alpha }}_{n,m}{\Vert }_2^2+\sum _{j=1}^{3p^2}{\lambda }_j|{\mathbf{\alpha }}_{n,m,j}|)\\ & \text{s.t.}\mathbf{D}=[{\mathbf{D}}_{\text{E}}{\mathbf{D}}_{\text{I}}],{\mathbf{D}}^{\top }\mathbf{D}=\mathbf{I},\end{array}$$

$$\text{vspace}-2mm{\lambda }_j=\lambda /(\sqrt{{\mathbf{S}}_k(j)}+\epsilon ),$$

$$\text{vspace}-2mm\begin{array}{r}{\mathbf{\alpha }}_{n,m}^{(t+1)}:=\underset{{\mathbf{\alpha }}_{n,m}}{\mathrm{argmin}}\Vert {\overline{\mathbf{y}}}_{n,m}-{\mathbf{D}}^{(t)}{\mathbf{\alpha }}_{n,m}{\Vert }_2^2+\sum _{j=1}^{3p^2}{\lambda }_j|{\mathbf{\alpha }}_{n,m,j}|.\end{array}$$

$${\mathbf{\alpha }}_{n,m}^{(t+1)}=\text{sgn}(({\mathbf{D}}^{(t)}{)}^{\top }{\overline{\mathbf{y}}}_{n,m})\odot \text{max}(|({\mathbf{D}}^{(t)}{)}^{\top }{\overline{\mathbf{y}}}_{n,m}|-\mathbf{\lambda },0),$$

$$\begin{array}{rl}{\mathbf{D}}_{\text{I}}^{(t+1)}:& =\underset{{\text{D}}_{\text{I}}}{\mathrm{argmin}}\sum _{n=1}^N\sum _{m=1}^M\Vert {\overline{\mathbf{y}}}_{n,m}-\mathbf{D}{\mathbf{\alpha }}_{n,m}^{(t+1)}{\Vert }_2^2\\ & =\underset{{\text{D}}_{\text{I}}}{\mathrm{argmin}}\Vert {\overline{\mathbf{Y}}}_n-\mathbf{D}{\mathbf{A}}^{(t+1)}{\Vert }_F^2\\ \text{s.t.}\mathbf{D}& =[{\mathbf{D}}_{\text{E}}{\mathbf{D}}_{\text{I}}],{\mathbf{D}}_{\text{I}}^{\top }{\mathbf{D}}_{\text{I}}={\mathbf{I}}_{(3p^2-r)},{\mathbf{D}}_{\text{E}}^{\top }{\mathbf{D}}_{\text{I}}=0,\end{array}$$

$$\text{vspace}-1mm(\mathbf{I}-{\mathbf{D}}_{\text{E}}{\mathbf{D}}_{\text{E}}^{\top })\mathbf{Y}({\mathbf{A}}_{\text{I}}^{(t+1)}{)}^{\top }={\mathbf{U}}_{\text{I}}{\mathbf{S}}_{\text{I}}{\mathbf{V}}_{\text{I}}^{\top }.$$

$$\begin{array}{rl}\hat{\mathcal{D}}=& \arg \underset{\mathcal{D}}{\min }\Vert \mathcal{Y}-\mathcal{D}\mathcal{A}{\Vert }_F^2\\ & \text{s.t.}{\mathcal{D}}^{\top }\mathcal{D}={\mathbf{I}}_{(3p^2-r)\times (3p^2-r)},{\mathcal{E}}^{\top }\mathcal{D}=0_{r\times (3p^2-r)},\end{array}$$

$$({\mathbf{I}}_{3p^2\times 3p^2}-\mathcal{E}{\mathcal{E}}^{\top })\mathcal{Y}{\mathcal{A}}^{\top }=\mathcal{U}\Sigma {\mathcal{V}}^{\top }$$

$$\text{vspace}-2mm{\hat{\mathbf{y}}}_{n,m}={\mathbf{D}}^{(T)}{\hat{\mathbf{\alpha }}}_{n,m}^{(T)}+{\mathbf{\mu }}_n,$$

$$\text{vspace}-2mm{\min }_{\mathbf{\alpha }}\Vert \mathbf{y}-\mathbf{D}\mathbf{\alpha }{\Vert }_2^2+{\sum }_{j=1}^{3p^2}{\lambda }_j|{\mathbf{\alpha }}_j|.$$

$$\text{vspace}-2mm{\min }_{\mathbf{\alpha }}\Vert {\mathbf{D}}^T\mathbf{y}-\mathbf{\alpha }{\Vert }_2^2+{\sum }_{j=1}^{3p^2}{\lambda }_j|{\mathbf{\alpha }}_j|.$$

$$\text{vspace}-2mm{\min }_{\mathbf{\alpha }}{\sum }_{j=1}^{3p^2}(({\mathbf{z}}_j-{\mathbf{\alpha }}_j{)}^2+{\lambda }_j|{\mathbf{\alpha }}_j|).$$

$$\text{vspace}-2mm{\min }_{{\mathbf{\alpha }}_j}({\mathbf{z}}_j-{\mathbf{\alpha }}_j{)}^2+{\lambda }_j|{\mathbf{\alpha }}_j|,$$

$$\text{vspace}-1mm\hat{\mathbf{\alpha }}=\text{sgn}({\mathbf{D}}^{\mathbf{T}}\mathbf{y})\odot \text{max}(|{\mathbf{D}}^{\mathbf{T}}\mathbf{y}|-\mathbf{\lambda },0),$$

$$\begin{array}{rl}\hat{\mathcal{D}}=& \arg {\min }_{\mathcal{D}}\Vert \mathcal{Y}-\mathcal{D}\mathcal{A}{\Vert }_F^2\\ & \text{s.t.}{\mathcal{D}}^{\top }\mathcal{D}={\mathbf{I}}_{(3p^2-r)\times (3p^2-r)},{\mathcal{E}}^{\top }\mathcal{D}=0_{r\times (3p^2-r)}.\end{array}$$

$$\begin{array}{rl}\hat{\mathcal{D}}& =\arg {\min }_{\mathcal{D}}\Vert \mathcal{Y}-\mathcal{D}\mathcal{A}{\Vert }_F^2=\arg {\max }_{\mathcal{D}}\text{Tr}(\mathcal{A}{\mathcal{Y}}^{\top }\mathcal{D})\\ & \text{s.t.}{\mathcal{D}}^{\top }\mathcal{D}={\mathbf{I}}_{(3p^2-r)\times (3p^2-r)},{\mathcal{E}}^{\top }\mathcal{D}=0_{r\times (3p^2-r)}.\end{array}$$

$$\partial \mathcal{L}/\partial \mathcal{D}=\mathcal{Y}{\mathcal{A}}^{\top }-\mathcal{D}({\Gamma }_1+{\Gamma }_1^{\top })-\mathcal{E}{\Gamma }_2^{\top }=0_{3p^2\times (3p^2-r)}.$$

$${\mathcal{E}}^{\top }\mathcal{Y}{\mathcal{A}}^{\top }={\Gamma }_2^{\top }.$$

$$({\mathbf{I}}_{3p^2\times 3p^2}-\mathcal{E}{\mathcal{E}}^{\top })\mathcal{Y}{\mathcal{A}}^{\top }=\mathcal{D}({\Gamma }_1+{\Gamma }_1^{\top }).$$

$$({\mathbf{I}}_{3p^2\times 3p^2}-\mathcal{E}{\mathcal{E}}^{\top })\mathcal{Y}{\mathcal{A}}^{\top }{\mathcal{D}}^{\top }=\mathcal{D}({\Gamma }_1+{\Gamma }_1^{\top }){\mathcal{D}}^{\top }.$$

$$({\mathbf{I}}_{3p^2\times 3p^2}-\mathcal{E}{\mathcal{E}}^{\top })\mathcal{Y}{\mathcal{A}}^{\top }{\mathcal{D}}^{\top }=\mathcal{U}\Sigma {\mathcal{V}}^{\top }{\mathcal{D}}^{\top }=\mathcal{D}\mathcal{V}\Sigma {\mathcal{U}}^{\top }.$$

$$\begin{array}{r}\hat{\mathcal{D}}=\arg {\min }_{\mathcal{D}}\Vert \mathcal{Y}-\mathcal{D}\mathcal{A}{\Vert }_F^2=\arg {\max }_{\mathcal{D}}\text{Tr}({\mathcal{Y}}^{\top }\mathcal{D}\mathcal{A}).\end{array}$$

$$({\mathbf{I}}_{3p^2\times 3p^2}-\mathcal{E}{\mathcal{E}}^{\top })\mathcal{Y}{\mathcal{A}}^{\top }={\mathcal{U}}^{\ast }{\mathcal{K}}_u\Sigma {\mathcal{K}}_v({\mathcal{V}}^{\ast }{)}^{\top },$$

$$\mathcal{U}\Sigma {\mathcal{V}}^{\top }{\mathcal{D}}^{\top }=\mathcal{D}\mathcal{V}\Sigma {\mathcal{U}}^{\top }\text{vspace}-2mm$$

$$\Sigma ={\mathcal{U}}^{\top }\mathcal{D}\mathcal{V}\Sigma {\mathcal{U}}^{\top }\mathcal{D}\mathcal{V}\text{vspace}-2mm$$

$$\Delta_{ii}=\left\{\begin{array}[]{ll}1&\textrm{if $1\leq i\leq\text{rank}(\Sigma)$};\\ \pm 1&\textrm{if $\text{rank}(\Sigma)<i\leq 3p^{2}-r$}.\\ \end{array}\right.\vspace{-1mm}$$

Full text

I Introduction

Image denoising is a crucial and indispensable step to improve image quality in digital imaging systems. In particular, with the decrease of size of CMOS/CCD sensors, image is more easily to be corrupted by noise and hence denoising is becoming increasingly important for high resolution imaging. The problem of image denoising has been extensively studied in literature and numerous image denoising methods [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44] have been proposed in the past decades. Most of existing denoising methods focus on the scenario of additive white Gaussian noise (AWGN) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 17, 18, 19, 15, 16, 20, 21, 22, 23, 24, 25], where the observed noisy image $\mathbf{y}$ is modeled as the addition of clean image $\mathbf{x}$ and AWGN $\mathbf{n}$, i.e., $\mathbf{y}=\mathbf{x}+\mathbf{n}$. There are also methods proposed for removing Poisson noise [26, 27], mixed Poisson and Gaussian noise [28, 29, 30, 31], mixed Gaussian and impulse noise [32, 33, 34], and realistic noise in real photography [35, 36, 37, 38, 41, 39, 40, 42, 44, 43].

Natural images have many properties, such as sparsity and nonlocal self-similarity, which can be employed as useful priors for designing image denoising methods. Based on the facts that natural images will be sparsely distributed in some transformed domain, wavelet [1] and curvelet [2] transforms have been widely adopted for image denoising. The sparse representation based methods [3, 4, 5, 6, 7, 8] encode image patches over a dictionary by using $\ell_{1}$-norm minimization to enforce the sparsity. The well-known bilateral filters [9] employ the prior information that image pixels exhibit similarity in both spatial domain and intensity domain. Other image priors such as multiscale self-similarity [10] and nonlocal self-similarity [11, 12], or the combination of multiple image priors [13, 14] have also been successfully used in image denoising. For example, by using low-rank minimization to characterize the image nonlocal self-similarity, the WNNM [13] method achieves state-of-the-art performance for AWGN denoising.

Instead of using predefined image priors, methods have also been proposed to learn priors from natural images for denoising. The generative image prior learning methods usually learn prior models from a set of external clean images and apply the learned prior models to the given noisy image [15, 16, 14, 17, 18, 19], or learn priors from the given noisy image to perform denoising [3]. Recently, the discriminative image prior learning methods [20, 21, 22, 25, 23, 24], which learn denoising models from pairs of clean and noisy images, have been becoming popular. The representative methods include the neural network based methods [20, 21, 22], random fields based methods [23, 24], and reaction diffusion based methods [25].

Most of the above mentioned methods focus on AWGN removal, however, the assumption of AWGN is too ideal to be true for real-world noisy images, where the noise is much more complex and varies with different scenes, cameras and camera settings (ISO, shutter speed, and aperture, etc.) [42, 45]. As a result, many denoising methods in literature, including those learning based methods, become less effective when applied to real-world noisy images. Fig. 1 shows an example, where we apply some representative and state-of-the-art denoising methods, including CBM3D [7], WNNM [13], DnCNN [22], CSF [24], and TNRD [25] to a real-world noisy image (captured by a Nikon D800 camera with ISO is 3200) provided in [42]. One can see that these methods either remain much the noise or over-smooth the image details.

There have been a few methods [35, 36, 37, 38, 42, 41, 39, 40, 43] and software toolboxes [44] developed for real-world noisy image denoising. Almost all of these methods follow a two-stage framework: first estimate the parameters of the noise model (usually assumed to be Gaussian or mixture of Gaussians (MoG)), and then perform denoising with the estimated noise model. However, the noise in real-world noisy images is very complex and is hard to be modeled by explicit distributions such as Gaussian and MoG. According to [45], the noise corrupted in the in-camera imaging process [46, 47, 42, 48] is signal dependent and comes from five main sources: photon shot, fixed pattern, dark current, readout, and quantization noise. The existing methods [35, 36, 37, 38, 42, 41, 39, 40, 44] mentioned above may not perform well on real-world noisy image denoising tasks. Fig. 1 also shows the denoising results of two real-world noisy image denoising methods, Noise Clinic [39, 40] and Neat Image [44]. One can see that these two methods still generate much noise caused artifacts.

This work aims to develop a new paradigm for real-world noisy image denoising. Different from existing real-world noisy image denoising methods [35, 36, 37, 38, 42, 41, 39, 40] which focus on noise modeling, we focus on image prior learning. We argue that with a strong and adaptive prior learning scheme, robust denoising performance on real-world noisy images can still be obtained. To achieve this goal, we propose to first learn image priors from external clean images, and then employ the learned external priors to guide the learning of internal priors from the given noisy image. The flowchart of the proposed method is illustrated in Fig. 2. We first extract millions of patch groups (PGs) from a set of high quality natural images, with which a Gaussian Mixture Model (GMM) is learned as the external image prior. The learned GMM prior model is used to assign each PG extracted from the given noisy image into its most suitable cluster via maximum a-posterior, and then an external-internal hybrid orthogonal dictionary is learned as the final prior for each cluster, with which the denoising can be readily performed by weighted sparse coding with closed form solution. The external priors learned from clean images preserve fine-scale image structural information, which is hard to be reproduced from noisy images. Therefore, external dictionary can serve as a good supplement to the internal dictionary. Our proposed denoising method is simple and efficient, yet our extensive experiments on real-world noisy images demonstrate its better denoising performance than the current state-of-the-arts.

II Related Work

II-A Internal and External Prior Learning

Learning natural image priors plays a key role in image denoising [10, 8, 3, 4, 5, 15, 16, 14, 17, 18, 19, 20, 21, 22, 23, 24, 25]. There are mainly four categories of prior learning based methods. 1) External prior learning methods [15, 16, 14] learn priors (e.g., dictionaries) from a set of external clean images, and the learned priors are used to recover the latent clean image from the given noisy image. 2) Internal prior learning methods [10, 8, 3, 4, 5] directly learn priors from a given noisy image, and image denoising is often done simultaneously with the prior learning process. 3) Discriminative prior learning methods [20, 21, 22, 23, 24, 25] learn discriminative models or mapping functions from clean and noisy image pairs, and the learned models or mapping functions are applied to a noisy image for denoising. 4) Hybrid methods [17, 18, 19] combine the external and internal priors to denoise the given input image.

It has been shown [15, 16, 14] that the external priors learned from natural clean images are effective and efficient for universal image denoising problems, whereas they are not adaptive to the given noisy image and some fine-scale image structures may not be well recovered. By contrast, the internal priors learned from the given noisy image are adaptive to image content, but the learned priors can be much affected by noise and the learning processing is usually slow [10, 8, 3, 4, 5]. Besides, most of the internal prior learning methods [10, 8, 3, 4, 5] assume additive white Gaussian noise (AWGN), making the learned priors less robust for real-world noisy images. In this paper, we use external priors to guide the internal prior learning. Our method is not only much faster than the traditional internal learning methods, but also very robust to denoise real-world noisy images.

In [17], the authors employed external clean patches to denoise noisy patches with high individual Signal-to-Noise-Ratio (PatchSNR), and employed internal noisy patches to denoise noisy patches with low PatchSNR. This is essentially different from our work which employs the external patch group based prior to guide the clustering and dictionary learning of the internal noisy patch groups. In [18], the external priors are only used to guide the internal patch clustering for image denoising, while in our work, the learned external priors are employed to guide not only the internal clustering, but also the internal dictionary learning. Besides, the method of [18] follows a patch based framework for AWGN removal, while in our work we employ a patch group based framework for real-world noisy image denoising. In addition, some technical details are also different. For example, method in [18] utilizes low-rank minimization for denoising, while we use dictionary learning and sparse coding for denoising. In the Targeted Image Denoising (TID) method [19], targeted images are selected from a large dataset for each patch in the input noisy image for denoising, which is computationally expensive.

II-B Real-World Noisy Image Denoising

Most of the denoising methods in literature [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 20, 21, 22, 23, 24, 25] assume AWGN noise and use simulated noisy images for algorithm design and evaluation. Recently, several denoising methods have been proposed to remove unknown noise from real-world noisy images [35, 36, 37, 38, 42, 41, 39, 40]. Portilla [35] employed a correlated Gaussian model to estimate the noise of each wavelet subband. Rabie [36] modeled the noisy pixels as outliers and performed denoising via Lorentzian robust estimator. Liu et al. [37] proposed the “noise level function” to estimate the noise and performed denoising by learning a Gaussian conditional random field. Gong et al. [38] proposed to model the data fitting term via weighted sum of $\ell_{1}$ and $\ell_{2}$ norms and performed denoising by a simple sparsity regularization term in the wavelet transform domain. The “Noise Clinic” [39, 40] estimates the noise distribution by using a multivariate Gaussian model and removes the noise by using a generalized version of nonlocal Bayesian model [12]. Zhu et al. [41] proposed a Bayesian method to approximate and remove the noise via a low-rank mixture of Gaussians (MoG) model. The method in [42] models the cross-channel noise in real-world noisy image as a multivariate Gaussian and the noise is removed by the Bayesian nonlocal means filter [49]. The commercial software Neat Image [44] estimates the noise parameters from a flat region of the given noisy image and filters the noise correspondingly.

The methods [35, 36, 37, 38, 42, 41, 39, 40] emphasize much on the noise modeling, and they use Gaussian or MoG to model the noise in real-world noisy images. Nonetheless, the noise in real-world noisy images is very complex and hard to be modeled by explicit distributions [45]. These works ignore the importance of learning image priors, which actually can be easier to model compared with modeling the complex realistic noise. In this paper, we propose a simple yet effective image prior learning method for real-world noisy image denoising. Due to its strong prior modeling ability, the proposed method simply models the noise as locally Gaussian, and it achieves highly competitive performance on real-world noisy image denoising.

III External Prior Guided Internal Prior Learning for Image Denoising

In this section, we first describe the learning of external prior, and then describe in detail the guided internal prior learning method, followed by the denoising algorithm.

III-A Learn External Patch Group Priors

The nonlocal self-similarity based patch group (PG) prior learning [14] has proved to be very effective for image denosing. In this work, we extract PGs from natural clean images to learn external priors. A PG is a group of similar patches to a local patch. In our method, each local patch is extracted from a RGB image with patch size $p\times p\times 3$. We search the $M$ most similar (i.e., smallest Euclidean distance) patches to this local patch (including the local patch itself) in a $W\times W$ region around it. Each patch is stretched to a patch vector $\mathbf{x}_{m}\in\mathbb{R}^{3p^{2}\times 1}$ to form the PG, denoted by $\{\mathbf{x}_{m}\}_{m=1}^{M}$. The mean vector of this PG is $\bm{\mu}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{x}_{m}$, and the group mean subtracted PG is defined as $\mathbf{\overline{X}}\triangleq\{\mathbf{\overline{x}}_{m}=\mathbf{x}_{m}-\bm{\mu}\}_{m=1}^{M}$.

Assume that a number of $L$ PGs are extracted from a set of external natural images, and the $l$-th PG is $\mathbf{\overline{X}}_{l}\triangleq\{\mathbf{\overline{x}}_{l,m}\}_{m=1}^{M},l=1,...,L$. A Gaussian Mixture Model (GMM) is learned to model the PG prior. The overall log-likelihood function is

[TABLE]

The learning process is similar to the GMM learning in [14, 16, 50]. Finally, a GMM model with $K$ Gaussian components is learned, and the learned parameters include mixture weights $\{\pi_{k}\}_{k=1}^{K}$, mean vectors $\{\bm{\mu}_{k}\}_{k=1}^{K}$, and covariance matrices $\{\bm{\Sigma}_{k}\}_{k=1}^{K}$. Note that the mean vector of each cluster is naturally zero, i.e., $\bm{\mu}_{k}=\bm{0}$.

To better describe the subspace of each Gaussian component, we perform singular value decomposition (SVD) [51] on the covariance matrix:

[TABLE]

The eigenvector matrices $\{\bm{U}_{k}\}_{k=1}^{K}$ will be employed as the external orthogonal dictionary to guide the internal sub-dictionary learning in next sub-section. The singular values in $\bm{S}_{k}$ reflect the significance of the singular vectors in $\bm{U}_{k}$. They will also be utilized as prior weights for weighted sparse coding in our denoising algorithm.

III-B Guided Internal Prior Learning

After the external PG prior model is learned from external natural clean images, we employ it to guide the internal PG prior learning for a given real-world noisy image. The guidance lies in two aspects. First, the external prior will guide the subspace clustering [52, 53] of internal noisy PGs. Second, the external prior will guide the orthogonal dictionary learning of internal noisy PGs.

III-B1 Internal Subspace Clustering

Given a real-world noisy image $\mathbf{y}$, we extract $N$ (overlapped) local patches from it. Similar to the external prior learning stage, for the $n$-th ($n=1,...,N$) local patch we search its $M$ most similar (by Euclidean distance) patches around it to form a noisy PG, denoted by $\bm{Y}_{n}=\{\mathbf{y}_{n,1},...,\mathbf{y}_{n,M}\}$. Then the group mean of $\bm{Y}_{n}$, denoted by $\bm{\mu}_{n}$, is subtracted from each patch by $\bm{\overline{y}}_{n,m}\triangleq\mathbf{y}_{n,m}-\bm{\mu}_{n}$, leading to the mean subtracted noisy PG $\bm{\overline{Y}}_{n}\triangleq\{\bm{\overline{y}}_{n,m}\}_{m=1}^{M}$.

The external GMM prior models $\{\mathcal{N}(\bm{0},\bm{\Sigma}_{k})\}_{k=1}^{K}$ basically characterize the subspaces of natural high quality PGs. Therefore, we project each noisy PG $\bm{\overline{Y}}_{n}$ into the subspaces of $\{\mathcal{N}(\bm{0},\bm{\Sigma}_{k})\}_{k=1}^{K}$ and assign it to the most suitable subspace based on the posterior probability:

[TABLE]

for $k=1,...,K$. Then $\bm{\overline{Y}}_{n}$ is assigned to the subspace with the maximum a-posteriori (MAP) probability $\max_{k}P(k|\bm{\overline{Y}}_{n})$.

III-B2 Guided Orthogonal Dictionary Learning

Assume that we have assigned all the internal noisy PGs $\{\bm{\overline{Y}}_{n}\}_{n=1}^{N}$ to their corresponding most suitable subspaces in $\{\mathcal{N}(\bm{0},\bm{\Sigma}_{k})\}_{k=1}^{K}$. For the $k$-th subspace, the noisy PGs assigned to it are $\{\bm{\overline{Y}}_{k_{n}}\}_{n=1}^{N_{k}}$, where $\bm{\overline{Y}}_{k_{n}}=[\bm{\overline{y}}_{k_{n},1},...,\bm{\overline{y}}_{k_{n},M}]$ and $\sum_{k=1}^{K}N_{k}=N$. We propose to learn an orthogonal dictionary $\bm{D}_{k}$ from each set of PGs $\bm{\overline{Y}}_{k_{n}}$ to characterize the internal PG prior with the guidance of the corresponding external orthogonal dictionary $\bm{U}_{k}$ (Eq. (2)). The reasons that we learn orthogonal dictionaries are two-fold. Firstly, the PGs $\{\bm{\overline{Y}}_{k_{n}}\}_{n=1}^{N_{k}}$ are in a subspace of the whole space of all PGs; therefore, there is no necessary to learn a redundant over-complete dictionary to characterize it, while an orthonormal dictionary has naturally zero mutual incoherence [54]. Secondly, the orthogonality of dictionary can make the patch encoding in the testing stage very efficient, leading to an efficient denoising algorithm (please refer to sub-section III-C for more details).

We let the orthogonal dictionary $\bm{D}_{k}$ be

[TABLE]

where $\bm{D}_{k,\text{E}}=\bm{U}_{k}(:,1:r)\in\mathbb{R}^{3p^{2}\times r}$ is the external sub-dictionary and it includes the first $r$ most important eigenvectors of $\bm{U}_{k}$, and the internal sub-dictionary $\bm{D}_{k,\text{I}}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ is to be adaptively learned from the noisy PGs $\{\bm{\overline{Y}}_{k_{n}}\}_{n=1}^{N_{k}}$. The rationale to design $\bm{D}_{k}$ as a hybrid dictionary is as follows. The external sub-dictionary $\bm{D}_{k,\text{E}}$ is pre-trained from external clean data, and it represents the $k$-th latent subspace of natural images, which is helpful to reconstruct the common latent structures of images. However, $\bm{D}_{k,\text{E}}$ is general to all images but not adaptive to the given noisy image. Some fine-scale details specific to the given image may not be well characterized by $\bm{D}_{k,\text{E}}$. Therefore, we learn an internal sub-dictionary $\bm{D}_{k,\text{I}}$ to supplement $\bm{D}_{k,\text{E}}$. In other words, $\bm{D}_{k,\text{I}}$ is to reveal the latent subspace adaptive to the input noisy image, which cannot be effectively represented by $\bm{D}_{k,\text{E}}$.

For notation simplicity, in the following development we ignore the subspace index $k$ for $\bm{\overline{Y}}_{k_{n}}$ and $\bm{D}_{k}$, etc. The learning of hybrid orthogonal dictionary $\bm{D}$ is performed under the following weighted sparse coding framework:

[TABLE]

where $\bm{I}$ is the $3p^{2}$ dimensional identity matrix, $\bm{\alpha}_{n,m}$ is the sparse coding vector of the $m$-th patch $\bm{\overline{y}}_{n,m}$ in the $n$-th PG $\bm{\overline{Y}}_{n}$ and $\bm{\alpha}_{n,m,j}$ is the $j$-th element of $\bm{\alpha}_{n,m}$. $\lambda_{j}$ is the $j$-th regularization parameter defined as

[TABLE]

where $\bm{S}_{k}(j)$ is the $j$-th singular value of diagonal singular value matrix $\bm{S}_{k}$ (please refer to Eq. (2)) and $\varepsilon$ is a small positive number to avoid zero denominator. Note that $\bm{D}_{\text{E}}=\bm{U}_{k}$ if $r=3p^{2}$ and $\bm{D}_{\text{E}}=\varnothing$ if $r=0$.

In the dictionary learning model (5), we use the $\ell_{2}$ norm to model the representation residual of PGs. This is because the patches in those PGs have similar content, and we assume that the noise therein will have similar statistics, which can be roughly modeled as locally Gaussian. On the other hand, this will make the dictionary learning much easier to solve. We employ an alternating iterative approach to solve the optimization problem (5). Specifically, we initialize the orthogonal dictionary as $\bm{D}^{(0)}=\bm{U}_{k}$ and for $t=0,1,...,T-1$, and alternatively update $\bm{\alpha}_{n,m}$ and $\bm{D}_{\text{I}}$ as follows.

Updating Sparse Coding Coefficients: Given the orthogonal dictionary $\bm{D}^{(t)}$, we update each sparse coding vector $\bm{\alpha}_{n,m}$ by solving

[TABLE]

Since dictionary $\bm{D}^{(t)}$ is orthogonal, the problems (7) has a closed-form solution

[TABLE]

where $\bm{\lambda}=\frac{1}{2}[\lambda_{1},\lambda_{2},...,\lambda_{3p^{2}}]^{\top}$ is the vector of regularization parameter, $\text{sgn}(\bullet)$ is the sign function and $\odot$ means element-wise multiplication. The detailed derivation of Eq. (8) can be found in Appendix A.

Updating Internal Sub-dictionary: Given the sparse coding vectors $\{\bm{\alpha}_{n,m}^{(t+1)}\}$, we update the internal sub-dictionary by solving

[TABLE]

where $\textbf{A}^{(t+1)}=[\bm{\alpha}_{1,1}^{(t+1)},...,\bm{\alpha}_{1,M}^{(t)},...,\bm{\alpha}_{N,1}^{(t+1)},...,\bm{\alpha}_{N,M}^{(t+1)}]$ and $\bm{I}_{(3p^{2}-r)}$ is the $(3p^{2}-r)$ dimensional identity matrix. The sparse coefficients matrix can be written as $\bm{A}^{(t+1)}=[(\bm{A}_{\text{E}}^{(t+1)})^{\top}\ (\bm{A}_{\text{I}}^{(t+1)})^{\top}]^{\top}$ where the external part $\bm{A}_{\text{E}}^{(t+1)}\in\mathbb{R}^{r\times NM}$ and the internal part $\bm{A}_{\text{I}}^{(t+1)}\in\mathbb{R}^{(3p^{2}-r)\times NM}$ represent the coding coefficients of $\bm{Y}$ over external sub-dictionary $\bm{D}_{\text{E}}$ and internal sub-dictionary $\bm{D}_{\text{I}}^{(t)}$, respectively. According to the following Theorem 1, by setting $\mathcal{Y}=\bm{\overline{Y}}_{n}-\bm{D}_{\text{E}}\bm{A}_{\text{E}}^{(t+1)},\mathcal{E}=\bm{D}_{\text{E}},\mathcal{D}=\bm{D}_{\text{I}},\mathcal{A}=\bm{A}_{\text{I}}$, the problem (9) has a closed-form solution $\bm{D}_{\text{I}}^{(t+1)}=\bm{U}_{\text{I}}\bm{V}_{\text{I}}^{\top}$, where $\bm{U}_{\text{I}}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ and $\bm{V}_{\text{I}}\in\mathbb{R}^{(3p^{2}-r)\times(3p^{2}-r)}$ are the orthogonal matrices obtained by the following SVD [51]

[TABLE]

The orthogonality of internal sub-dictionary $\bm{D}_{\text{I}}^{(t+1)}$ can be checked by $(\bm{D}_{\text{I}}^{(t+1)})^{\top}(\bm{D}_{\text{I}}^{(t+1)})=\bm{V}_{\text{I}}\bm{U}_{\text{I}}^{\top}\bm{U}_{\text{I}}\bm{V}_{\text{I}}^{\top}=\bm{I}_{(3p^{2}-r)}$. In fact, the Theorem 1 provides a sufficient and necessary condition to guarantee the existence of the closed-form solution for the internal sub-dictionary of the problem (9).

Theorem 1**.**

Let $\mathcal{A}\in\mathbb{R}^{(3p^{2}-r)\times M}$, $\mathcal{Y}\in\mathbb{R}^{3p^{2}\times M}$ be two given data matrices. $\mathcal{E}\in\mathbb{R}^{3p^{2}\times r}$ is a given matrix satisfying $\mathcal{E}^{\top}\mathcal{E}=\bm{I}_{r\times r}$, then $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$ is the necessary condition of

[TABLE]

where $\mathcal{U}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ and $\mathcal{V}\in\mathbb{R}^{(3p^{2}-r)\times(3p^{2}-r)}$ are the orthogonal matrices obtained by performing economy (a.k.a. reduced) SVD [51]:

[TABLE]

Besides, if $\text{rank}(\Sigma)=3p^{2}-r$, $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$ is also the sufficient condition of problem (11).

The proof of the Theorem 1 can be found in Appendix B. Though the problem (9) has a closed-form solution by SVD [51], the uniqueness of solution cannot be guaranteed since the matrices $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}$ as well as $\mathcal{U}$ and $\mathcal{V}$ may be reduced to matrices of lower rank. Hence, we also analyze the uniqueness of the solution $\hat{\mathcal{D}}$ by the following Theorem 2, whose proof can be found in Appendix C.

Theorem 2**.**

(a) If $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ is nonsingular, i.e., $\text{rank}(\Sigma)=3p^{2}-r$, then the solution of $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$ is unique; (b) If $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}$ is singular, i.e., $0\leq\text{rank}(\bm{\Sigma})<3p^{2}-r$, then the number of possible solutions of $\hat{\mathcal{D}}$ is $2^{3p^{2}-r-\text{rank}(\bm{\Sigma})}$ for fixed $\mathcal{U}$ and $\mathcal{V}$.

The above alternative updating steps are repeated until the number of iterations exceeds a preset threshold. In each step, the energy value of the objective function (5) is decreased and we empirically found that the proposed model usually converges in 10 iterations. We summarize the procedures in Algorithm 1.

III-C The Denoising Algorithm

The denoising of the given noisy image $\mathbf{y}$ can be simultaneously done with the guided internal sub-dictionary learning process. Once we obtain the solutions of sparse coding vectors $\{\hat{\bm{\alpha}}_{n,m}^{(T)}\}$ in Eq. (8) and the orthogonal dictionary $\bm{D}^{(T)}=[\bm{D}_{\text{E}}\ \bm{D}_{\text{I}}^{(T)}]$ in Eq. (9), the latent clean patch $\hat{\mathbf{y}}_{n,m}$ of the $m$-th noisy patch in PG $\bm{Y}_{n}$ is reconstructed as

[TABLE]

where $\bm{\mu}_{n}$ is the group mean of $\bm{Y}_{n}$. The latent clean image is then reconstructed by aggregating all the reconstructed patches in all PGs. We perform the above denoising procedures for several iterations for better denoising outputs. The proposed denoising algorithm is summarized in Algorithm 2.

IV Experiments

IV-A Implementation Details

The noise in real-world images is very complex due to the many factors such as sensors, lighting conditions and camera settings. It is difficult to evaluate one algorithm by tuning its parameters for all these different settings. In this work, we fix the parameters of our algorithm and apply it to all the testing datasets, though they were captured by different types of sensors and under different camera settings. The parameters of our method include the patch size $p$, the number of similar patches $M$ in a patch group (PG), the window size $W$ for PG searching, the number of Gaussian components $K$ in GMM, the number of atoms $r$ in the external sub-dictionaries, the sparse regularization parameter $\lambda$, the iteration numbers $T$ for solving problem (5) and $IteNum$ for Alg. 2.

The performance of our proposed method varies little when we set patch size between $p=6$ and $p=9$, and we fix the patch size as $p=6$ to save computational cost. The search window is fixed to $W=31$ to balance computational cost and denoising accuracy of the proposed method. The number of patches in a patch group is set as $M=10$, while using more patches will not bring clear benefits. We learn the external GMM prior with 3.6 million PGs extracted from the Kodak PhotoCD Dataset (http://r0k.us/graphics/kodak/), which includes 24 high quality color images. The number of Gaussians in GMM is set as $K=32$, while using more Gaussians can only bring slightly better performance but cost more computational resources. The number of atoms in the external sub-dictionaries affects little the performance when it is set between $r=27$ and $r=81$, and we set it as $r=54$ to make the external and internal sub-dictionaries have the same number of atoms. We set the number of iterations as $T=2$ for solving the problem (5), while the number of iterations for Alg. 2 is set as $IteNum=4$.

One key parameter of our model is the regularization parameter $\lambda$. Fig. 3 plots the curves of PSNR/SSIM results w.r.t $\lambda$ on the 15 cropped image in dataset [42]. One can see that our proposed method achieves good PSNR/SSIM performance within a certain range of $\lambda$. Similar observations can be made on other datasets. We fix $\lambda=0.001$ in the paper, and it works well across the three datasets used in our experiments.

All the parameters of our method are fixed in all experiments, which are run under the Matlab2014b environment on a machine with Intel(R) Core(TM) i7-5930K CPU of 3.5GHz and 32GB RAM. We will release the code with the publication of this work.

IV-B The Testing Datasets

We evaluate the proposed method on three real-world noisy image datasets, where the images were captured under indoor or outdoor lighting conditions by different types of cameras and camera settings.

Dataset 1. The first dataset is provided in [42], which includes noisy images of 11 static scenes. The noisy images were collected under controlled indoor environment. Each scene was shot 500 times under the same camera and camera setting. The mean image of the 500 shots is roughly taken as the “ground truth”, with which the PSNR and SSIM [55]can be computed.

Since the image size is very large (about $7000\times 5000$) and the 11 scenes share repetitive contents, the authors of [42] cropped 15 smaller images (of size $512\times 512$) to perform experiments. In order to evaluate the proposed methods more comprehensively, we cropped 60 images of size $500\times 500$ from the dataset for experiments. Some samples are shown in Fig. 4. Note that our cropped 60 images and the 15 cropped images by the authors of [42] are from different shots.

Dataset 2 is called the Darmstadt Noise Dataset (DND) [56], which includes 50 different pairs of images of the same scenes captured by Sony A7R, Olympus E-M10, Sony RX100 IV, and Huawei Nexus 6P. The real-world noisy images are collected under higher ISO values with shorter exposure time, while the “ground truth” images are captured under lower ISO values with longer exposure times. Since the captured images are of megapixel-size, the authors cropped 20 bounding boxes of $512\times 512$ pixels from each image in the dataset, yielding $50\times 20=1000$ test crops in total. Some samples are shown in Fig. 5. Note that the “ground truth” images of this dataset have not been released yet, but one can submit the denoised images to the Project Website and get the average PSNR (dB) and SSIM results.

Dataset 3. On one hand, the scenes of Dataset 1 are mostly printed photos, and they cannot represent realistic objects and scenes with different reflectance properties. On the other hand, the Dataset 2 contains repetitive contents in the 20 cropped images for each of the 50 scenes. To remedy the limitations of Dataset 1 and Dataset 2, we construct another dataset which contains images of 10 different scenes captured by Canon 80D and Sony A7II cameras with more ISO settings and more comprehensive scenes. The ISO settings in our dataset are 800, 1600, 3200, 6400, 12800 while those of Dataset 1 are 1600, 3200, 6400. Compared to Dataset 2, our new dataset is more comprehensive on scene contents. Similar to Dataset 1, each scene was captured 500 shots, and the mean image of these 500 shots can be used a kind of ground-truth to evaluate the denoising algorithms. Fig. 6 shows some cropped images of the scenes in our dataset. One can see that the images contain a lot of different realistic objects with varying colors, shapes, materials, etc.

Our dataset provides real-world noisy images of realistic objects with different ISO settings. It can be used to more fairly evaluate the performance of different real-world noisy image denoising methods. Consider that the image resolution is very high (about $4000\times 4000$), for the convenience of experimental studies, we cropped 100 (10 for each scene) smaller images (of size $512\times 512$) from it to perform experiments. The whole dataset will be made publically available with the publication of this paper.

IV-C Comparison among external, internal and guided internal priors

To demonstrate the advantages of external prior guided internal prior learning, we perform real-world noisy image denoising by using external priors only (denoted by “External”), internal priors only (denoted by “Internal”), and the proposed guided internal priors (denoted by “Guided Internal”), respectively. For the “External” method, we utilize the full external dictionaries (i.e., $r=108$ in Eq. (5)) for denoising. For the “Internal” method, the overall framework is similar to the method of [5]. A GMM model (with $K=32$ Gaussians) is directly learned from the PGs extracted from the given noisy image without using any external data, and then the internal orthogonal dictionaries are obtained via Eq. (2) to perform denoising. All parameters of the “External” and “Internal” methods are tuned to achieve their best performance.

We compare the three methods on the 60 cropped images from Dataset 1 [42]. The average PSNR and run time are listed in Table I. The best results are highlighted in bold. It can be seen that “Guided Internal” method achieves better PSNR than both “External” and “Internal” methods. In addition, the “Internal” method is very slow because it involves online GMM learning, while the “Guided Internal” method is only a little slower than the “External” method. Figs. 7 and 8 show the denoised images of two noisy images by the three methods. One can see that the “External” method is good at recovering large-scale structures (see Fig. 7) while the “Internal” method is good at recovering fine-scale textures (see Fig. 8). By utilizing external priors to guide the internal prior learning, our proposed method can effectively recover both the large-scale structures and fine-scale textures.

IV-D Comparison with State-of-the-Art Denoising Methods

Comparison methods. We compare the proposed method with state-of-the-art image denoising methods, including GAT-BM3D [30], CBM3D [7], WNNM [13], TID [19], MLP [20], DnCNN [22], CSF [24], TNRD [25], Noise Clinic (NC) [39, 40], Cross-Channel (CC) [42], and Neat Image (NI) [44]. Among these methods, GAT-BM3D [30] is a state-of-the-art Poisson noise reduction method. The method CBM3D [7] is a state-of-the-art method for color image denoising and the noise on color images is assumed to be additive white Gaussian. The methods of WNNM, MLP, DnCNN, CSF, and TNRD are state-of-the-art Gaussian noise removal methods for grayscale images, and we apply them to each channel of color images for denoising. NC is a blind image denoising method, and NI is a set of commercial software for image denoising, which has been embedded into Photoshop and Corel PaintShop. The code of CC is not released but its results on the 15 cropped images are available at [42]. Therefore, we only compare with it on the 15 cropped images in Dataset 1 [42].

Noise level of comparison methods. For the CBM3D method, the standard deviation of noise on color images should be given as a parameter. For methods of WNNM, MLP, CSF, and TNRD, the noise level in each color channel should be input. For the DnCNN method, it is trained to deal with noise in a range of levels $0\sim 55$. We retrain the models of discriminative denoising methods MLP, CSF, and TNRD (using the released codes by the authors) at different noise levels from $\sigma=5$ to $\sigma=50$ with a gap of $5$. The denoising is performed by processing each channel with the model trained at the same (or nearest) noise level. The noise levels ($\sigma_{r},\sigma_{g},\sigma_{b}$) in R, G ,B channels are assumed to be Gaussian and can be estimated via some noise estimation methods [57, 58]. In this paper, we employ the method [58] to estimate the noise level for each color channel.

Results on Dataset 1. As described in section 4.2, there is a mean image for each of the 11 scenes used in Dataset 1 [42], and those mean images can be roughly taken as “ground truth” images for quantitative evaluation of denoising algorithms. We firstly perform quantitative comparison on the 15 cropped images used in [42]. The results on PSNR (dB) and speed (second) of GAT-BM3D, CBM3D, WNNM, TID, MLP, CSF, TNRD, DnCNN, NC, NI and CC are listed in Table II (The results of CC are copied from the original paper [42]). The best PSNR results of each image are highlighted in bold. One can see that on 8 out of the 15 images, our method achieves the best PSNR values. CC achieves the best PSNR on 3 of the 15 images. It should be noted that in the CC method, a specific model is trained for each camera and camera setting, while our method uses the same model for all images. On average, our proposed method has 0.27dB PSNR improvements over the second best method CC and much higher PSNR gains over other competing methods. The method GAT-BM3D does not work well on most images. This is because real world noise is much more complex than Poisson.

Figs. 9 and 10 show the denoised images of one scene captured by Canon 5D Mark 3 at ISO = 3200 and Nikon D800 at ISO = 6400, respectively. We can see that GAT-BM3D, CBM3D, TID, DnCNN, NC, NI and CC would either remain noise or generate artifacts, while TNRD over-smooths much the image. By using the external prior guided internal priors, our proposed method preserves edges and textures better than other methods while removing the noise, leading to visually more pleasant outputs. Specifically, Fig. 10 is used to illustrate the denoising performance of our method on fine-scale textures such as hair, which is a very challenging task. Even the “ground truth” mean image cannot show very clear details of the hair. Though our method cannot reproduce clearly the details (e.g., the local direction of hair in some regions), it demonstrates the best visual results among the competing methods. More comparisons on visual quality and SSIM [55] index can be found in the supplementary file.

We then perform denoising experiments on the 60 images we cropped from [42]. The average PSNR results are listed in Table III (CC is not compared since the code is not available). Again, our proposed method achieves much better PSNR results than the other methods. The improvements of our method over the second best method (TNRD) are 0.43dB on PSNR. Fig. 11 shows the denoised images of one scene captured by Nikon D800 at ISO = 3200. We can see again that the proposed method obtain better visual quality than other competing methods. More comparisons on visual quality and SSIM can be found in the supplementary file.

Results on Dataset 2. In Table IV, we list the average PSNR (dB) results of the competing methods on the 1000 cropped images in the DND dataset [56]. We can see again that the proposed method achieves better performance than the other competing methods. Note that the “ground truth” images of this dataset have not been released yet, so we are not able to calculate the PSNR and SSIM results for each noisy image in this dataset, nor compare with the “ground truth” mean image. However, one can submit the denoised images to the project website and get the average PSNR and SSIM results on the whole 1000 images. Fig. 12 shows the denoised images of a scene “0001_2” captured by a Nexus 6P phone [56]. The noise level in this image is relatively high. Hence, this image can be used to justify the performance of the proposed method on real-world noisy images with lower PSNR (around 20dB). One can see that the proposed method achieves visually more pleasing results than the other denoising methods. More comparisons on visual quality and SSIM can be found in the supplementary file.

Results on Dataset 3. Similar to Dataset 1 [42], there is a “ground truth” image for each of the 10 scenes used in our constructed Dataset 3. We perform quantitative comparison on the 100 cropped images. The average PSNR results of competing methods are listed in Table IV. We can see that our proposed method achieves much better PSNR results than the other methods. The improvements of our method over the second best method (TNRD) is 0.16dB on PSNR. Fig. 13 shows the denoised images of one scene captured by Canon 80D at ISO = 12800. We can see again that the proposed method removes the noise while maintains better details (such as the vertical black shadow area) than other competing methods. More comparisons on visual quality and SSIM can be found in the supplementary file.

Comparison on speed. Efficiency is an important aspect to evaluate the efficiency of algorithms. We compare the speed of all competing methods except for CC. All experiments are run under the Matlab2014b environment on a machine with Intel(R) Core(TM) i7-5930K CPU of 3.5GHz and 32GB RAM. The average running time (second) of the compared methods on the 100 real-world noisy images is shown in Table V. The least average running time are highlighted in bold. One can easily see that the commercial software Neat Image (NI) is the fastest method with highly optimized code. For a $512\times 512$ image, NI costs about 0.6 second. The other methods cost from 5.2 (TNRD) to 152.2 (WNNM) seconds, while the proposed method costs about 24.1 seconds. It should be noted that GAT-BM3D, CBM3D, TNRD, and NC are implemented with compiled C++ mex-function and with parallelization, while WNNM, TID, MLP, CSF, DnCNN, and the proposed method are implemented purely in Matlab.

V Conclusion

We proposed a new prior learning method for the real-world noisy image denoising problem by exploiting the useful information in both external and internal data. We first learned Gaussian Mixture Models (GMMs) from a set of clean external images as general image prior, and then employed the learned GMM model to guide the learning of adaptive internal prior from the given noisy image. Finally, a set of orthogonal dictionaries were output as the external-internal hybrid prior models for image denoising. Extensive experiments on three real-world noisy image datasets, including a new dataset constructed by us by different types of cameras and camera settings, demonstrated that our proposed method achieves much better performance than state-of-the-art image denoising methods in terms of both quantitative measure and visual perceptual quality.

Appendix A Closed-Form Solution of the Weighted Sparse Coding Problem (7)

For notation simplicity, we ignore the indices $n,m,t$ in problem (7). It turns into the following weighted sparse coding problem:

[TABLE]

Since $\bm{D}$ is an orthogonal matrix, problem (14) is equivalent to:

[TABLE]

For simplicity, we denote $\mathbf{z}=\mathbf{D^{T}y}$. Here we have $\lambda_{j}>0$, $j=1,...,3p^{2}$, then problem (15) can be written as:

[TABLE]

The problem (16) is separable w.r.t. each $\bm{\alpha}_{j}$ and hence can be simplified to $3p^{2}$ independent scalar minimization problems:

[TABLE]

where $j=1,...,3p^{2}$. Taking derivative of $\bm{\alpha}_{j}$ in problem (17) and setting the derivative to be zero. There are two cases for the solution.

(a) If $\bm{\alpha}_{j}\geq 0$, we have $2(\bm{\alpha}_{j}-\mathbf{z}_{j})+\lambda_{j}=0,$ and the solution is $\hat{\bm{\alpha}}_{j}=\mathbf{z}_{j}-\frac{\lambda_{j}}{2}\geq 0.$ So $\mathbf{z}_{j}\geq\frac{\lambda_{j}}{2}>0$, and the solution $\hat{\bm{\alpha}}_{j}$ can be written as $\hat{\bm{\alpha}}_{j}=\text{sgn}(\mathbf{z}_{j})*(|\mathbf{z}_{j}|-\frac{\lambda_{j}}{2}),$ where $\text{sgn}(\bullet)$ is the sign function.

(b) If $\bm{\alpha}_{j}<0$, we have $2(\bm{\alpha}_{j}-\mathbf{z}_{j})-\lambda_{j}=0$ and the solution is $\hat{\bm{\alpha}}_{j}=\mathbf{z}_{j}+\frac{\lambda_{j}}{2}<0.$ So $\mathbf{z}_{j}<-\frac{\lambda_{j}}{2}<0$, and the solution $\hat{\bm{\alpha}}_{j}$ can be written as $\hat{\bm{\alpha}}_{j}=\text{sgn}(\mathbf{z}_{j})*(-\mathbf{z}_{j}-\frac{\lambda_{j}}{2})=\text{sgn}(\mathbf{z}_{j})*(|\mathbf{z}_{j}|-\frac{\lambda_{j}}{2}).$

In summary, we have the final solution of the weighted sparse coding problem (14) as:

[TABLE]

where $\bm{\lambda}=\frac{1}{2}[\lambda_{1},\lambda_{2},...,\lambda_{3p^{2}}]^{\top}$ is the vector of regularization parameter and $\odot$ means element-wise multiplication.

Appendix B Proof of the Theorem 1

Let $\mathcal{A}\in\mathbb{R}^{(3p^{2}-r)\times M},\mathcal{Y}\in\mathbb{R}^{3p^{2}\times M}$ be two given data matrices. Denote by $\mathcal{E}\in\mathbb{R}^{3p^{2}\times r}$ the external subdictionary and $\mathcal{D}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ the internal subdictionary. For simplicity, we assume $3p^{2}\geq M$. The problem in Theorem 1 is as follows:

[TABLE]

Proof.

We firstly prove the necessary condition. Since $\mathcal{D}^{\top}\mathcal{D}=\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}$, we have

[TABLE]

The Lagrange function is $\mathcal{L}=\text{Tr}(\mathcal{A}\mathcal{Y}^{\top}\mathcal{D})-\text{Tr}(\Gamma_{1}(\mathcal{D}^{\top}\mathcal{D}-\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}))-\text{Tr}(\Gamma_{2}(\mathcal{D}^{\top}\mathcal{E}))$, where $\Gamma_{1}$ and $\Gamma_{2}$ are the Lagrange multipliers. Take the derivative of $\mathcal{L}$ w.r.t. $\mathcal{D}$ and set it to be matrix $\bm{0}$ of conformal dimensions, we can get

[TABLE]

Since $\mathcal{D}^{\top}\mathcal{D}=\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}$ and $\mathcal{E}^{\top}\mathcal{D}=\bm{0}_{3p^{2}\times(3p^{2}-r)}$, by left multiplying both sides of the Eq. (22) by $\mathcal{E}^{\top}$, we have

[TABLE]

Put the Eq. (22) back into Eq. (21), we have

[TABLE]

Right multiplying both sides of Eq. (23) by $\mathcal{D}^{\top}$, we have

[TABLE]

This shows that $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}\mathcal{D}^{\top}$ is a symmetric matrix of order $3p^{2}\times 3p^{2}$. Then we perform economy (or reduced) singular value decomposition (SVD) [51] on $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}=\mathcal{U}\Sigma\mathcal{V}^{\top}$, there is

[TABLE]

Hence, we have $\mathcal{U}=\mathcal{D}\mathcal{V}$, or equivalently $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$. The necessary condition is proved.

Now we prove the sufficient condition. If $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$, then $\hat{\mathcal{D}}^{\top}\hat{\mathcal{D}}=\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}$. To prove $\mathcal{E}^{\top}\hat{\mathcal{D}}=\bm{0}_{3p^{2}\times(3p^{2}-r)}$, we left multiply both sides of Eq. (25) by $\mathcal{E}^{\top}$ and have $\bm{0}_{3p^{2}\times(3p^{2}-r)}=\mathcal{E}^{\top}(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}\hat{\mathcal{D}}^{\top}=\mathcal{E}^{\top}\mathcal{U}\Sigma\mathcal{V}^{\top}\hat{\mathcal{D}}^{\top}=\mathcal{E}^{\top}\mathcal{U}\Sigma\mathcal{U}^{\top}$ . It means that $\mathcal{E}^{\top}\mathcal{U}\Sigma\mathcal{U}^{\top}=\bm{0}_{3p^{2}\times 3p^{2}}$. This only happens when $\mathcal{E}^{\top}\mathcal{U}=\bm{0}_{3p^{2}\times(3p^{2}-r)}$ since $\text{rank}(\Sigma)=3p^{2}-r$ and $\mathcal{U}\Sigma\mathcal{U}^{\top}$ is positive definite. Then $\mathcal{E}^{\top}\hat{\mathcal{D}}=\mathcal{E}^{\top}\mathcal{U}\mathcal{V}^{\top}=\bm{0}_{3p^{2}\times(3p^{2}-r)}$.

Finally we prove that $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$ is the solution of

[TABLE]

Note that by cyclic perturbation which retains the trace unchanged and due to $\mathcal{E}^{\top}\hat{\mathcal{D}}=\bm{0}_{3p^{2}\times(3p^{2}-r)}$, we have $\text{Tr}(\mathcal{Y}^{\top}\hat{\mathcal{D}}\mathcal{A})=\text{Tr}(\mathcal{Y}\mathcal{A}^{\top}\hat{\mathcal{D}}^{\top})=\text{Tr}((\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}\hat{\mathcal{D}}^{\top})=\text{Tr}(\mathcal{U}\Sigma\mathcal{V}^{\top}\mathcal{V}\mathcal{U}^{\top})=\text{Tr}(\Sigma).$ For every $\mathcal{D}$ satisfying that $\mathcal{D}^{\top}\mathcal{D}=\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}$, $\mathcal{E}^{\top}\mathcal{D}=\bm{0}_{3p^{2}\times(3p^{2}-r)}$, we have $\text{Tr}(\mathcal{Y}^{\top}\mathcal{D}\mathcal{A})=\text{Tr}((\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}\mathcal{D}^{\top})=\text{Tr}(\mathcal{U}\Sigma\mathcal{V}^{\top}\mathcal{D}^{\top})=\text{Tr}(\Sigma\mathcal{V}^{\top}\mathcal{D}^{\top}\mathcal{U})$. By using a generalization version [59] of the Kristof’s Theorem [60], we have $\text{Tr}(\mathcal{Y}^{\top}\mathcal{D}\mathcal{A})=\text{Tr}(\Sigma\mathcal{V}^{\top}\mathcal{D}^{\top}\mathcal{U})\leq\text{Tr}(\Sigma).$ The equality is obtained at $\mathcal{V}^{\top}\mathcal{D}^{\top}\mathcal{U}=\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}$, i.e., $\mathcal{D}=\mathcal{U}\mathcal{V}^{\top}=\hat{\mathcal{D}}$. This completes the proof. ∎

Appendix C Proof of the Theorem 2

Before we prove the Theorem 2, we need firstly prove the following Lemma 1.

Lemma 1: Let $\mathcal{E}\in\mathbb{R}^{3p^{2}\times r}$ be an orthogonal matrix with $\mathcal{E}^{\top}\mathcal{E}=\bm{I}_{r\times r}$, then $\text{rank}(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\geq 3p^{2}-r$.

Proof.

Since $\text{rank}(\mathcal{E}\mathcal{E}^{\top})\leq\min\{\text{rank}(\mathcal{E}),\text{rank}(\mathcal{E}^{\top})\}=r$ and $\text{rank}(\mathcal{E}\mathcal{E}^{\top})\geq\text{rank}(\mathcal{E})+\text{rank}(\mathcal{E}^{\top})-r=2r-r=r$ by Sylvester’s inequality, we have $\text{rank}(\mathcal{E}\mathcal{E}^{\top})=r$. Then, $\text{rank}(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\geq\text{rank}(\bm{I}_{3p^{2}\times 3p^{2}})-\text{rank}(\mathcal{E}\mathcal{E}^{\top})\geq 3p^{2}-r$. ∎

The $\text{rank}(\Sigma)$ ($\Sigma$ is defined in Theorem 1) depends on $\text{rank}(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})$, $\text{rank}(\mathcal{Y})$ and $\text{rank}(\mathcal{A})$. Note that $\text{rank}(\mathcal{Y})\geq M$ and $\text{rank}(\mathcal{A})\geq\min\{3p^{2},M\}$ and $\text{rank}(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\geq 3p^{2}-r$. Hence, $\text{rank}(\Sigma)\leq\min\{3p^{2}-r,M\}$.

Now we prove the Theorem 2:

Proof.

a) If $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ is nonsingular, i.e., $\text{rank}(\Sigma)=3p^{2}-r$, $\Sigma$ may have distinct or multiple non-zero singular values. In the SVD [51] of $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}=\mathcal{U}\Sigma\mathcal{V}^{\top}$, the singular vectors in $\mathcal{U}$ and $\mathcal{V}$ can be determined up to orientation. Hence, we can reformulate it as

[TABLE]

where $\mathcal{U}^{*}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$ and $\mathcal{V}^{*}\in\mathbb{R}^{(3p^{2}-r)\times(3p^{2}-r)}$ are arbitrarily orientated singular vectors of $\mathcal{U}$ and $\mathcal{V}$, respectively. The $\mathcal{K}_{u}$ and $\mathcal{K}_{v}$ are diagonal matrices with $+1$ or $-1$ as diagonal elements in arbitrary distribution. $\Sigma\in\mathbb{R}^{(3p^{2}-r)\times(3p^{2}-r)}$ is a diagonal matrix with singular values in non-increasing order, i.e., $\Sigma_{11}\geq\Sigma_{22}\geq...\geq\Sigma_{(3p^{2}-r)(3p^{2}-r)}\geq 0$. If we fix $\mathcal{K}_{u}$, then $\mathcal{K}_{v}$ is uniquely determined to meet the above requirements of $\Sigma$. If the orientations of the singular vectors of $\mathcal{U}^{*}$ are fixed, then $\mathcal{U}=\mathcal{U}^{*}\mathcal{K}_{u}$ is determined, so do the orientations of the singular vectors of $\mathcal{V}^{*}$ and $\mathcal{V}^{\top}=\mathcal{K}_{v}(\mathcal{V}^{*})^{\top}$. In this case, the solution of $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}=\mathcal{U}^{*}\mathcal{K}_{u}\mathcal{K}_{v}(\mathcal{V}^{*})^{\top}$ is unique. When $\Sigma$ has multiple singular values, the unique solution of $\hat{\mathcal{D}}$ can be proved in a similar way.

b) If $(\bm{I}_{3p^{2}\times 3p^{2}}-\mathcal{E}\mathcal{E}^{\top})\mathcal{Y}\mathcal{A}^{\top}$ is singular, i.e., $0\leq\text{rank}(\Sigma)<3p^{2}-r$, and $\Sigma$ has $3p^{2}-r-\text{rank}(\Sigma)$ (at least one) zero singular values. The discussion in a) can still be applied to the singular vectors corresponding to the nonzero singular values, and the production of these singular vectors in $\mathcal{U}$ and $\mathcal{V}$ is still unique. However, the singular vectors corresponding to the zero singular values could be in arbitrary orientations as long as they satisfy the conditions of $\mathcal{U}^{\top}\mathcal{U}=\mathcal{V}^{\top}\mathcal{V}=\mathcal{V}\mathcal{V}^{\top}=\bm{I}_{(3p^{2}-r)\times(3p^{2}-r)}$. Since $\mathcal{U}\in\mathbb{R}^{3p^{2}\times(3p^{2}-r)}$, $\mathcal{U}\mathcal{U}^{\top}$ no longer equals to the identity matrix of order $3p^{2}\times 3p^{2}$. From Eq. (25), we have

[TABLE]

Right multiplying both sides of Eq. (28) by $\mathcal{D}\mathcal{V}$ and left multiplying each side by $\mathcal{U}^{\top}$, we have

[TABLE]

Hence, $\Delta=\mathcal{U}^{\top}\mathcal{D}\mathcal{V}\in\mathbb{R}^{(3p^{2}-r)\times(3p^{2}-r)}$ is a diagonal matrix, the diagonal elements of which are

[TABLE]

Thus, we have $\mathcal{D}=\mathcal{U}\Delta\mathcal{V}^{\top}$. That is, if $\text{rank}(\Sigma)<3p^{2}-r$, once we get the solution of $\hat{\mathcal{D}}=\mathcal{U}\mathcal{V}^{\top}$ in problem (19), $\mathcal{D}=\mathcal{U}\Delta\mathcal{V}^{\top}$ with suitable $\Delta$ is also the solution of problem (19). In fact, the number of solutions $\hat{\mathcal{D}}$ for problem (19) is $2^{3p^{2}-r-\text{rank}(\Sigma)}$ given fixed $\mathcal{U}$ and $\mathcal{V}$. ∎