Color Image and Multispectral Image Denoising Using Block Diagonal Representation

TL;DR

This paper introduces a novel patch-level block diagonal representation for color and multispectral image denoising, demonstrating that a trained global basis combined with local PCA transforms yields competitive results efficiently.

Contribution

It emphasizes the importance of patch-level representation using block diagonal matrices and proposes a simple, effective transform-threshold-inverse method with fast implementation.

Findings

Competitive denoising results on simulated and real datasets

Robustness and efficiency demonstrated through extensive experiments

Effective use of global patch basis and local PCA transforms

Abstract

Filtering images of more than one channel is challenging in terms of both efficiency and effectiveness. By grouping similar patches to utilize the self-similarity and sparse linear approximation of natural images, recent nonlocal and transform-domain methods have been widely used in color and multispectral image (MSI) denoising. Many related methods focus on the modeling of group level correlation to enhance sparsity, which often resorts to a recursive strategy with a large number of similar patches. The importance of the patch level representation is understated. In this paper, we mainly investigate the influence and potential of representation at patch level by considering a general formulation with block diagonal matrix. We further show that by training a proper global patch basis, along with a local principal component analysis transform in the grouping dimension, a simple transform-threshold-inverse method could produce very competitive results. Fast implementation is also developed to reduce computational complexity. Extensive experiments on both simulated and real datasets demonstrate its robustness, effectiveness and efficiency.

Figures

Tables

Table 1. TABLE I: Techniques and priors adopted by different methods based on the grouping-collaborative filtering-aggregation framework. LR: Low Rank; SC: Sparse Coding; HT: Hard-thresholding; WF: Wiener Filter [ 2 ]

Method	LLRT [21]	WTR1 [51]	TWSC [15]	NCSR [52]	CBM3D	Proposed
Techique	LR	LR	SC	SC	HT+WF	HT

Table 2. TABLE II: Brief description of three real-world color image datasets.

Dataset	Type	# of Images	Image Size
CC15 [28]	Real-World	15	$512\times 512$
CC60 [32]	Real-World	60	$500\times 500$
Xu-100 [43]	Real-World	100	$512\times 512$
Ours [46]	Real-World	249	$1024\times 1024$

Table 3. TABLE III: Input noise level of Gaussian denoisers for CC15, CC60 and Xu-100 datasets.

Method	LSCD	WTR1	FFD-Net	LLRT
$\sigma$	10	15	15	20
Method	4DHOSVD1	CBM3D	MSt-SVD	CMSt-SVD
$\sigma$	30	20	30	25

Table 4. TABLE IV: PSNR results and average computational time (s) of different methods on CC15 dataset. The three best results are bolded.

Methods	Canon 5D ISO = 3200			Canon D600 ISO = 3200			Nikon D800 ISO = 1600			Nikon D800 ISO = 3200			Nikon D800 ISO = 6400			Average	Time (s)
LSCD	37.86	36.21	35.52	34.65	36.26	38.24	37.90	38.88	38.32	37.45	36.49	37.73	32.33	32.55	32.62	36.20	9.68
LLRT	39.23	36.31	35.93	34.74	36.83	40.58	37.39	40.27	37.78	39.79	37.34	41.03	35.09	34.05	34.11	37.36	$>1000$
WTR1	41.09	36.92	36.25	34.68	36.48	40.52	38.26	41.40	38.61	39.98	37.70	41.36	35.16	34.22	34.43	37.81	$>2000$
TNRD	39.51	36.47	36.45	34.79	36.37	39.49	38.11	40.52	38.17	37.69	35.90	38.21	32.81	32.33	32.29	36.61	N/A
MLP	39.00	36.34	36.33	34.70	36.20	39.33	37.95	40.23	37.94	37.55	35.91	38.15	32.69	32.33	32.29	36.46	N/A
DnCNN	37.26	34.13	34.09	33.62	34.48	35.41	35.79	36.08	35.48	34.08	33.70	33.31	29.83	30.55	30.09	33.86	N/A
GID	40.82	37.19	36.92	35.32	36.62	38.68	38.88	40.66	39.20	37.92	36.62	37.64	33.01	32.93	32.96	37.02	55.60
WCMNNM	41.20	37.25	36.48	35.54	37.03	39.56	39.26	41.45	39.54	38.94	37.40	39.42	34.85	33.97	33.97	37.72	318.29
FFD-Net	39.40	37.02	36.53	34.97	36.73	41.02	38.66	41.53	38.80	40.15	37.61	41.18	34.13	33.66	33.69	37.68	28.98
TWSC	40.55	35.92	35.15	35.36	37.09	41.13	39.36	41.91	38.81	40.27	37.22	42.09	35.53	34.15	33.93	37.90	480.80
CBM3D	40.77	37.31	36.98	35.21	36.76	40.13	39.02	41.65	39.40	39.59	37.49	39.47	34.13	33.73	33.85	37.69	6.98
CBM3D_best	40.96	37.31	37.15	35.38	36.81	40.45	39.25	41.65	39.59	39.86	37.54	40.38	34.85	33.92	34.16	37.95	6.98
4DHOSVD1	40.22	36.97	36.55	35.02	36.60	39.78	38.85	41.35	39.11	39.24	37.28	39.47	34.40	33.81	34.01	37.51	120.18
MSt-SVD	40.33	37.25	36.83	35.16	36.71	40.29	38.97	41.49	39.24	39.61	37.43	39.93	34.34	33.82	33.96	37.69	110.06
CMSt-SVD	40.79	37.37	37.01	35.29	36.95	40.93	39.21	41.98	39.54	39.98	37.65	40.05	34.50	33.93	34.01	37.95	98.88

Table 5. TABLE V: Average PSNR and SSIM values of compared methods on CC60 and Xu-100 datasets. The best results are bolded.

Dataset	Index	LLRT	WTR1	FFD-Net	GID	TWSC	MCWNNM	4DHOSVD1	CBM3D	CBM3D_best	CMSt-SVD
CC60	PSNR	38.51	39.69	39.73	38.41	39.66	39.03	39.15	39.40	39.68	39.75
	SSIM	0.9636	0.9764	0.9770	0.9633	0.9759	0.9698	0.9729	0.9740	0.9775	0.9756
Xu	PSNR	38.51	38.56	38.56	38.37	38.62	38.51	38.51	38.69	38.81	38.82
	SSIM	0.9707	0.9669	0.9658	0.9675	0.9674	0.9671	0.9673	0.9694	0.9712	0.9694

Table 6. TABLE VI: Average PSNR and SSIM values of compared methods on our dataset. The best results are bolded. The input noise level of compared methods are tuned to produce their best average results of every camera.

Camera	# of Images	Index	LLRT	FFD-Net	GID	TWSC	MCWNNM	4DHOSVD1	CBM3D	CBM3D_best	CMSt-SVD
HUAWEI HONOR 6X	30	PSNR	39.54	40.05	39.52	39.71	39.46	39.82	39.97	40.48	40.08
		SSIM	0.9669	0.9669	0.9653	0.9651	0.9610	0.9658	0.9669	0.9740	0.9674
IPHONE 5S	36	PSNR	40.02	40.60	40.12	40.27	39.87	40.68	40.77	41.25	40.84
		SSIM	0.9676	0.9645	0.9642	0.9617	0.9567	0.9664	0.9668	0.9758	0.9668
IPHONE 6S	67	PSNR	39.72	40.49	40.16	40.12	40.18	40.36	40.55	41.16	40.53
		SSIM	0.9663	0.9707	0.9670	0.9619	0.9628	0.9671	0.9693	0.9783	0.9674
CANON 100D	55	PSNR	41.84	41.67	40.86	41.65	41.47	41.41	41.69	42.08	41.99
		SSIM	0.9784	0.9768	0.9743	0.9767	0.9774	0.9771	0.9780	0.9808	0.9794
CANON 600D	25	PSNR	42.53	42.55	41.60	42.52	42.07	42.14	42.54	42.89	42.75
		SSIM	0.9816	0.9824	0.9790	0.9824	0.9795	0.9810	0.9836	0.9851	0.9840
SONY A6500	36	PSNR	45.71	45.71	44.94	45.48	45.37	45.56	45.70	45.81	45.89
		SSIM	0.9899	0.9901	0.9887	0.9896	0.9894	0.9901	0.9902	0.9904	0.9903

Table 7. TABLE VII: Quantitative results comparison under i.i.d Gaussian noise σ = { 10 , 30 , 50 , 100 } 𝜎 10 30 50 100 \sigma=\{10,30,50,100\} on CAVE dataset.

Methods	$\sigma =10$				$\sigma =30$				$\sigma =50$				$\sigma =100$
	PSNR	SSIM	ERGAS	SAM	PSNR	SSIM	ERGAS	SAM	PSNR	SSIM	ERGAS	SAM	PSNR	SSIM	ERGAS	SAM
Noisy	28.13	0.4371	236.40	0.7199	18.59	0.1069	676.01	1.0085	14.15	0.0475	1126.7	1.1461	8.13	0.0136	2253.4	1.3074
PARAFAC	35.39	0.8759	108.78	0.2363	33.65	0.8294	125.13	0.3291	31.52	0.7393	154.19	0.4351	27.13	0.4637	250.68	0.6681
LRTA	41.36	0.9499	49.53	0.1718	36.06	0.8775	90.50	0.2446	33.52	0.8201	121.15	0.2897	30.06	0.7138	180.03	0.3649
LRMR	39.27	0.9094	64.81	0.3343	31.36	0.6451	157.65	0.6021	26.67	0.4000	264.28	0.7534	26.67	0.1850	469.26	0.9306
4DHOSVD1	45.43	0.9811	30.83	0.1084	39.78	0.9336	59.12	0.2272	36.82	0.8722	83.36	0.3385	32.66	0.7307	134.34	0.5599
BM4D	44.61	0.9784	33.32	0.1289	38.80	0.9283	65.23	0.2579	35.98	0.8685	91.19	0.3557	31.84	0.7197	144.91	0.5160
TDL	44.41	0.9797	34.32	0.1048	39.07	0.9493	63.18	0.1493	36.46	0.9171	85.24	0.2008	32.92	0.8284	128.15	0.3132
ISTReg	45.77	0.9802	30.53	0.1086	40.51	0.9488	53.05	0.1374	37.75	0.9271	70.16	0.1619	33.01	0.8648	120.77	0.2376
LLRT	46.60	0.9868	26.75	0.0842	41.49	0.9681	48.50	0.1221	38.65	0.9482	67.56	0.1551	35.39	0.9154	99.37	0.1962
MSt-SVD	45.20	0.9814	32.05	0.1064	40.23	0.9530	56.26	0.1737	37.73	0.9285	75.03	0.2231	34.20	0.8800	115.24	0.3142

Table 8. TABLE VIII: Quantitative results comparison under non i.i.d Gaussian noise with input noise level σ = 36 𝜎 36 \sigma=36 for all compared methods on CAVE dataset.

Index	PARAFAC	LRTA	4DHOSVD1	BM4D	TDL	LLRT	MSt-SVD
PSNR	32.7	33.0	38.1	36.6	32.3	38.3	39.3
SSIM	0.79	0.74	0.91	0.86	0.74	0.92	0.94
ERGAS	137.9	132.6	71.9	89.7	140.9	70.8	62.4
SAM	0.40	0.45	0.33	0.40	0.48	0.21	0.19
Time (m)	4.8	0.6	5.5	10.1	1.2	48.8	6.2

Table 9. TABLE IX: Choice of parameters for 4DHOSVD and MSt-SVD on color image and multispectral image (MSI) experiments.

Parameter	4DHOSVD		MSt-SVD
	Color Image	MSI	Color Image	MSI
$ps$	8	8	8	8
$K$	30	30	30	30
$SR$	20	16	20	16
$N_{step}$	4	4	4	4
$\gamma$	0.8	1	1.1 (), 1.2 ($\sigma >30$)	4

Equations

$$\hat{\mathbf{G}}=\arg \underset{{\mathbf{G}}_c}{\min }\Vert {\mathbf{G}}_n-{\mathbf{G}}_c{\Vert }_F^2+\rho \cdot \Psi ({\mathbf{G}}_c)$$

$$\hat{\mathbf{C}}=\arg \underset{C}{\min }\Vert {\mathbf{G}}_n-\mathbf{DC}{\Vert }_F^2+\lambda \Vert \mathbf{C}{\Vert }_1$$

$$\mathcal{T}({\mathcal{G}}_{ht})=\{\begin{array}{r}\mathcal{T}({\mathcal{G}}_n),|\mathcal{T}({\mathcal{G}}_n)|\ge \tau \\ 0,|\mathcal{T}({\mathcal{G}}_n)|<\tau \end{array}$$

$${\hat{p}}_i=\sum _{{\hat{p}}_{i_k}\in \hat{\mathbf{G}}}{w}_{i_k}{\hat{p}}_{i_k}$$

$$\mathcal{C}=\mathcal{G}{\times }_1{\mathbf{U}}_{row}^T{\times }_2{\mathbf{U}}_{col}^T{\times }_3{\mathbf{U}}_{color}^T{\times }_4{\mathbf{U}}_{group}^T$$

$$\mathbf{U}_{color}^{T}=\left(\begin{array}[]{ccc}\textbf{1/3}&\textbf{1/3}&\textbf{1/3}\\ 0.5&0&-0.5\\ 0.25&-0.5&0.25\\ \end{array}\right)$$

$$\begin{array}{rl}& \min \Vert \mathcal{G}-\mathcal{C}{\times }_1{\mathbf{U}}_{row}{\times }_2{\mathbf{U}}_{col}{\times }_3{\mathbf{U}}_{color}{\times }_4{\mathbf{U}}_{group}{\Vert }_F^2\\ & s.t{\mathbf{U}}_{row}^T{\mathbf{U}}_{row}=\mathbf{I},{\mathbf{U}}_{col}^T{\mathbf{U}}_{col}=\mathbf{I}\\ & {\mathbf{U}}_{color}^T{\mathbf{U}}_{color}=\mathbf{I},{\mathbf{U}}_{group}^T{\mathbf{U}}_{group}=\mathbf{I}\end{array}$$

$$\mathcal{C}=(\mathcal{G}{\times }_3{\mathbf{U}}_{color}^T){\times }_1{\mathbf{U}}_{row}^T{\times }_2{\mathbf{U}}_{col}^T{\times }_4{\mathbf{U}}_{group}^T$$

$$(\mathcal{G}{\times }_3{\mathbf{U}}_{color}^T{)}_i={\mathcal{P}}_i{\times }_3{\mathbf{U}}_{color}^T,i=1,2,...,K$$

$$bdiag(\hat{\mathcal{P}}_{i})=\left(\begin{array}[]{ccc}\hat{\mathcal{P}}_{i}(:,:,1)&&\\ &\hat{\mathcal{P}}_{i}(:,:,2)&\\ &&\hat{\mathcal{P}}_{i}(:,:,3)\\ \end{array}\right)$$

$${\hat{\mathcal{P}}}_i(:,:,k)=\sum _{j=1}^3{\mathbf{U}}_{color}(j,k){\mathcal{P}}_i(:,:,j),k=1,2,3$$

$$fdiag(\hat{\mathcal{G}})=\left(\begin{array}[]{ccc}\hat{\mathcal{G}}_{i}(:,:,1,:)&&\\ &\hat{\mathcal{G}}_{i}(:,:,2,:)&\\ &&\hat{\mathcal{G}}_{i}(:,:,3,:)\\ \end{array}\right)$$

$$\begin{array}{rl}& \mathcal{C}=fdiag(\hat{\mathcal{G}}){\times }_{1}bdiag({\mathbf{U}}_{row}^T){\times }_{2}bdiag({\mathbf{U}}_{col}^T){\times }_3{\mathbf{U}}_{group}^T\end{array}$$

$$bdiag(\mathbf{U})=\left(\begin{array}[]{ccc}\mathbf{U}&&\\ &\mathbf{U}&\\ &&\mathbf{U}\\ \end{array}\right)$$

$$bcirc(\mathcal{P}_{i})=\left(\begin{array}[]{ccc}\mathbf{P}_{i}(:,:,1)&\mathbf{P}_{i}(:,:,3)&\mathbf{P}_{i}(:,:,2)\\ \mathbf{P}_{i}(:,:,2)&\mathbf{P}_{i}(:,:,1)&\mathbf{P}_{i}(:,:,3)\\ \mathbf{P}_{i}(:,:,3)&\mathbf{P}_{i}(:,:,2)&\mathbf{P}_{i}(:,:,1)\\ \end{array}\right)$$

$$bcirc(\mathcal{G})=\left(\begin{array}[]{ccc}\mathcal{G}(:,:,1,:)&\mathcal{G}(:,:,3,:)&\mathcal{G}(:,:,2,:)\\ \mathcal{G}(:,:,2,:)&\mathcal{G}(:,:,1,:)&\mathcal{G}(:,:,3,:)\\ \mathcal{G}(:,:,3,:)&\mathcal{G}(:,:,2,:)&\mathcal{G}(:,:,1,:)\\ \end{array}\right)$$

$${\mathcal{C}}_{bcirc}=bcirc(\mathcal{G}){\times }_1{\mathbf{U}}_{bcir{c}_{row}}^T{\times }_2{\mathbf{U}}_{bcir{c}_{col}}^T{\times }_3{\mathbf{U}}_{group}^T$$

$${\mathcal{C}}_{fdiag}=fdiag(\hat{\mathcal{G}}){\times }_1{\mathbf{U}}_{fdia{g}_{row}}^T{\times }_2{\mathbf{U}}_{fdia{g}_{col}}^T{\times }_3{\mathbf{U}}_{group}^T$$

$$\left(\begin{array}[]{ccc}\textbf{1}&\textbf{1}&\textbf{1}\\ 1&-0.5-0.8660i&-0.5+0.8660i\\ 1&-0.5+0.8660i&-0.5-0.8660i\\ \end{array}\right)$$

$$bdiag(\hat{\mathcal{C}})=bdiag(\hat{\mathcal{A}})bdiag(\hat{\mathcal{B}})$$

$$\mathcal{A}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^T$$

$$bdiag(\hat{\mathcal{A}})=bdiag(\hat{\mathcal{U}})bdiag(\hat{\mathcal{S}})bdiag(\hat{\mathcal{V}}{)}^T$$

$$\mathcal{C}=\mathcal{G}{\ast }_1{\mathcal{U}}_{row}^T{\ast }_2{\mathcal{U}}_{column}^T{\ast }_3{\mathcal{U}}_{group}^T$$

$$\begin{array}{rl}& \min \sum _{i=1}^K\Vert {\mathcal{P}}_i-{\mathcal{U}}_{row}\ast {\mathcal{S}}_i\ast {\mathcal{U}}_{column}^T\Vert \\ & s.t{\mathcal{U}}_{row}^T\ast {\mathcal{U}}_{row}=\mathcal{I},{\mathcal{U}}_{column}^T\ast {\mathcal{U}}_{column}=\mathcal{I}\end{array}$$

$${\mathcal{G}}_{filtered}=\mathcal{C}{\ast }_1{\mathcal{U}}_{row}{\ast }_2{\mathcal{U}}_{column}{\ast }_3{\mathcal{U}}_{group}$$

$$\Vert {\mathcal{P}}_i-{\mathcal{P}}_j{\Vert }_F=\sqrt{3}\Vert bcirc({\mathcal{P}}_i)-bcirc({\mathcal{P}}_j){\Vert }_F$$

$${\mathbf{U}}_{group}={\mathbf{U}}_{bcirc\_group}$$

$$bdiag(\hat{\mathcal{P}})=(\mathbf{F}\otimes \mathbf{I})bcirc(\mathcal{P})(\mathbf{F}\otimes \mathbf{I}{)}^{-1}$$

$$\Vert {\mathcal{C}}_{fdiag}{\Vert }_F=\Vert {\mathcal{C}}_{bcirc}{\Vert }_F$$

Full text

I Introduction

IMAGE denoising plays an important role in modern image processing systems. The past few decades witness great achievements in this field [1], and methods that utilize self-similarity and non-local characteristics of natural images have drawn much attention due to their simplicity and effectiveness. Recently, great achievement is made by the well-known BM3D algorithm [2] which combines the nonlocal filters [3] and transform domain techniques [4]. Methods that share similar idea and procedure of BM3D are widely adopted to handle grayscale image [5, 6, 7]. When the input is color image (sRGB) or multispectral image that contains rich information and delivers more faithful representation for real scenes, directly applying the grayscale denoising algorithm to each channel [3] or spectral band [8] often fails to produce satisfactory results, and therefore efforts to understand and address noise reduction issue have been made from several different perspectives.

First, two main solutions are proposed to improve the channel-by-channel or band-by-band approach. The first strategy proposes to transform the original image into a less correlated color or band space, such that denoising in each transformed channel or band could be performed independently. the representative work is the color BM3D (CBM3D) which applies BM3D to the luminance-chrominance space [9]. The second strategy jointly characterizes the RGB channels or bands for better use of spectral correlation. Methods that fall into such category are widely considered with different priors and regularization. Briefly, [10] proposes a multichannel nonlocal fusion (MNLF) approach, [11] introduces color line to model the correlation among neighbouring pixels and channels, and [12] considers the spatial and spectral dependencies. Sparse and low rank priors are also adopted in several competitive methods [13, 14, 15]. Besides, to avoid vectorization of image patch, some recent works incorporate tensor representation [16] with higher order singular value decomposition (HOSVD) [17, 18, 7], low rank tensor approximation (LRTA) framework [19], Laplacian Scale Mixture modeling [20] and Hyper-Laplacian regularization [21].

In addition to the design of denoising strategy, noise modeling is also important. Most of existing methods consider additive white Gaussian noise (AWGN) and some efficient noise estimation methods [22, 23] can be employed. Besides, some non i.i.d. Gaussian denoisers are proposed for filtering Poission noise [24], mixed Gaussian and impulsive noise [25] and stripe noise [26]. In fact, noise in real-world images may be multiplicative and signal dependent [27], making noise modeling and estimation much more complex and challenging. [28] and [29] consider the non-linear processing steps in the camera pipeline in the noise model, and [30, 31, 32] combine external and internal priors. Some methods [33, 34, 35, 36], including the well-known software toolbox Neat Image (NI)111Neatlab ABSoft. https://ni.neatvideo.com/home are developed for noise reduction of real-world images. Apart from the conventional transform-domain approaches, many recent competitive methods [37, 38, 39, 40] are based on the advent of deep learning technique as a powerful feature extraction tool.

In order to compare different denoising methods, several real-world color image and multispectral image datasets [41, 28, 42, 43, 44, 45, 46] of various scenes are constructed, and each scene of a color image includes noisy and ”ground-truth” image pairs. A simple and reasonable approach adopted by [41, 28, 43, 46] to obtain ”ground-truth” image is to capture the same and unchanged scene for many times and compute their mean image. Different from the image averaging approach, [42] utilizes the Tobit regression to estimate the parameters of the noise process by accessing only two images, and [44] generates a high quality smartphone image denoising dataset with careful post-processing. Interestingly, recent experiments on real-world datasets show that BM3D based methods [9, 47] still demonstrate the most competitive performance in terms of both effectiveness and efficiency. The implementation details are unknown and there is debate [48, 49] that they may have touched the ceiling of image denoising.

Many competitive methods [13, 21, 15] attempt to approach the optimal performance by modeling the redundancy and correlation at group level with some iterative strategies [50] and a large number of similar patches. However, influence of the patch level representation is less carefully studied. Although the use of tensor representation may help preserve some structure information, the straightforward folding and unfolding operation may not fully exploit the relationship among all channels or spectral bands. In this paper, we investigate the potential and influence of patch level representation, and establish a general formulation with block diagonal matrix. We demonstrate that the combination of a proper global patch basis and local PCA can produce very competitive performance in terms of both efficiency and effectiveness. Extension to non-Gaussian noise in hyperspectral image is also discussed. Efforts have been made to reduce computational complexity, and all results reported in our paper could be reproduced very efficiently.

The rest of this paper is organized as follows: In Section II, related work is studied. In Section III, general formulation and the proposed denoising method are introduced. Section IV presents experiments on both simulated and real-world datasets. Conclusions are drawn in Section V.

II Related Works and Formulation

II-A Notations

Tensor is a multidimensional array, also known as a multi-way array, and its order is defined as the number of its dimension. In this paper, we mainly adopt the mathematical notations and preliminaries of tensors from [16]. Vectors and matrices are first- and second- order tensors which are denoted by boldface lowercase letters $\mathbf{a}$ and capital letters $\mathbf{A}$, respectively. A higher order tensor (the tensor of order three or above) is denoted by calligraphic letters, e.g., $\mathcal{A}$. An $N$th-order tensor is denoted as $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$. The $n$-mode product of a tensor $\mathcal{A}$ by a matrix $\mathbf{U}\in\mathrm{R}^{P_{n}\times I_{n}}$, denoted by $\mathcal{A}\times_{n}\mathbf{U}$ is also a tensor. The mode-$n$ matricization or unfolding of $\mathcal{A}$, denoted by $\mathbf{A}_{(n)}$, maps tensor elements $(i_{1},i_{2},\ldots,i_{N})$ to matrix element $(i_{n},j)$ where $j=1+\sum_{k=1,k\neq n}^{N}(i_{k}-1)J_{k}$, with $J_{k}=\prod_{m=1,m\neq n}^{k-1}I_{m}$. The Frobenius norm of a tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ is defined as $\|\mathcal{A}\|_{F}=\sqrt{\sum_{i_{1}=1}...\sum_{i_{N}=1}\mathcal{A}_{i_{1}...i_{N}}^{2}}$.

II-B Framework and Problem Formulation

The most popular and successful framework credited to [2] follows three consecutive steps: grouping, collaborative filtering and averaging. The flowchart of this paradigm is illustrated as Fig. 1. Specifically, given a reference patch $\mathcal{P}_{ref}$, the grouping step stacks some similar overlapping patches located in a local window $\Omega_{SR}$ into a group represented by matrix $\mathbf{G}$ or higher order tensor $\mathcal{G}$ with certain matching criteria [53, 54, 55]. One simple and commonly adopted metric is Euclidean distance measured by $\|\mathcal{P}_{ref}-\mathcal{P}\|_{F}^{2}$, $\forall\mathcal{P}\in\Omega_{SR}$. Collaborative filtering is then performed on group $\mathbf{G}$ to utilize the nonlocal similarity feature and estimate clean underlying patches from noisy observation, and it can be generally formulated as

[TABLE]

where $\mathbf{G}_{n}$ and $\mathbf{G}_{c}$ are noisy and underlying clean group of patches, respectively, $\|\mathbf{G}_{n}-\mathbf{G}_{c}\|_{F}^{2}$ measures the conformity between $\mathbf{G}_{c}$ and $\mathbf{G}_{n}$, and $\Psi(\mathbf{G}_{c})$ represents certain priors [56]. Low rank approximation is adopted in [57, 13, 21, 51] based on nuclear norm minimization [58] with $\Psi(\mathbf{G}_{c})=\|\mathbf{G}_{c}\|_{\ast}$, or tensor trace norm [59] with $\Psi(\mathcal{G}_{c})=\sum_{n=1}^{N}\alpha_{n}\|\mathbf{G_{c}}_{(n)}\|_{\ast}$. Authors in [60, 15, 52, 61] utilize sparse coding scheme that represents $\hat{\mathbf{G}}$ with a dictionary $\mathbf{D}$ and sparse coding atoms $\mathbf{C}$ by minimizing

[TABLE]

The state-of-the-art BM3D and HOSVD algorithms attempt to model sparsity in the transform domain by shrinking coefficients $\mathcal{T}(\mathcal{G}_{n})$ under a pre-defined threshold $\tau$ via

[TABLE]

Some representative techniques and priors are listed in Table I. After collaborative filtering, the estimated clean patches are averagely written back to their original location to further smooth out noise. More specifically, every pixel $\hat{p}_{i}$ of the denoised image is the (weighted) average of all pixels at the same position of filtered group $\hat{\mathbf{G}}$, which can be formulated as

[TABLE]

where $w_{i_{k}}$ and $\hat{p}_{i_{k}}$ denote weight and pixel, respectively.

II-C Multiway filtering technique

In this paper, we focus on the collaborative filtering step with transform domain technique, and regardless of different categorization criteria, it can be roughly modeled as an inverse problem with multiway filtering technique [62], which applies a specific transform and constraint to each dimension of a group of similar image patches. The organization of image patches and choice of transform and constraint often determine the denoising strategy, and the variations are quite extensive [9, 7, 20, 51]. In this subsection we introduce the representative 4DHOSVD and CBM3D for color image denoising, and our analysis can be extended to multispectral image with slight modification. For simplicity, a group of similar patches $\mathcal{P}$ is denoted by $\mathcal{G}$.

Both CBM3D and 4DHOSVD can be represented with the fourth-order tensor decomposition framework in [16] via

[TABLE]

where $\mathcal{C}$ is core tensor (coefficient), $\mathbf{U}_{row}$, $\mathbf{U}_{col}$, $\mathbf{U}_{color}$, and $\mathbf{U}_{group}$ are corresponding mode transforms. The major difference between CBM3D and 4DHOSVD is that CBM3D uses pre-defined discrete cosine transform (DCT) and opponent color space represented as a $3\times 3$ matrix

[TABLE]

specifically, the first slice of each patch in the new color space can be regarded as luminance channel, and the other two slices as chrominance channel. CBM3D is very efficient because it does not have to train local transforms, and its grouping step is performed only on the luminance channel. For 4DHOSVD, all mode transform matrices including $\mathbf{U}_{color}$ are learned by solving

[TABLE]

Compared with CBM3D, 4DHOSVD imposes a stronger constraint and requires the orthogonality of color space transform.

III Image Denoising Using Block Diagonal Representation

Using a 4D transform in equation (5) for CBM3D may be a little confusing, because after a certain color space transform, the original R, G, B channels are computed separatedly in the new color space, which also holds for 4DHOSVD if all mode transforms are obtained. Therefore, it may be expected to re-formulate (5) as independent channel-wise (third-order tensor) transform. In this section, we first generalize patch-level representation via block diagonal matrix, then discuss the proper choice for patch-level basis, and explain how it could be properly incorporated into the block diagonal representation and efficiently applied to image denoising.

III-A Block Diagonal Representation and Formulation of f-diagonal Tensor Decomposition

We notice that (5) can be rewritten as

[TABLE]

where $\mathcal{G}\times_{3}\mathbf{U}_{color}^{T}$ is equivalent to applying $\mathbf{U}_{color}$ to each patch of $\mathcal{G}$ via

[TABLE]

where $K$ and $\mathcal{P}_{i}$ represent the number of similar patches and the $i$-th patch of $\mathcal{G}$, respectively. For simplicity, we use $\hat{\mathcal{P}}_{i}$ to denote $\mathcal{P}_{i}\times_{3}\mathbf{U}_{color}^{T}$. Then we could define the block diagonal operator $bdiag(\hat{\mathcal{P}}_{i})$ via

[TABLE]

where each matrix on the diagonal position is a linear combination of all frontal slices of $\mathcal{P}_{i}$ via

[TABLE]

according to (9) and (10), $\mathcal{G}\times_{3}\mathbf{U}_{color}^{T}$ can be denoted as an f-diagonal tensor $fdiag(\hat{\mathcal{G}})$

[TABLE]

where $\hat{\mathcal{G}}=\mathcal{G}\times_{3}\mathbf{U}_{color}^{T}$. Based on (8) and (12), the 4D transform (5) is equivalent to

[TABLE]

where we make an abuse use of (10) to denote $bdiag(\mathbf{U}_{row})$ and $bdiag(\mathbf{U}_{col})$ as

[TABLE]

The same group representation $\mathbf{U}_{group}$ is applied for all frontal slices in (12) mainly because of two reasons. First, the patch-wise similarity often used in the grouping process does not guarantee slice-wise simialrity. Furthermore, according to equation (11), the new color space are a linear combination of R, G, B channels, thus are not totally de-correlated. Therefore, there exists a trade-off between the slice-wise and patch-wise relationship. Interestingly, CBM3D subtly takes care of this issue by considering only the luminance channel similarity. To narrow such gap, a suitable alternative to utilizing more group-level information (grouping more patches), is the recursive use of patch-level correlation via block circulant representation (BCR) [63].

Specifically, for each patch $\mathcal{P}_{i}\in\mathbb{R}^{ps\times ps\times 3}$ of color image, its BCR $bcirc(\mathcal{P}_{i})$ is a block circulant matrix [63] of size $3ps\times 3ps$ defined by

[TABLE]

thus, similar to (12), a block circulant tensor $bcirc(\mathcal{G})$ of size $3ps\times 3ps\times K$ can be denoted as

[TABLE]

Fig. 2 gives a straightforward illustration of (15) and (16), and some interesting feature is given in Appendix A.

Following the idea of multiway transform technique in (5) and based on (27), the third-order tensor decomposition of (16) is

[TABLE]

Obtaining two factor matrices $\mathbf{U}_{bcirc_{row}}$ and $\mathbf{U}_{bcirc_{col}}$ requires the eigenvalue decomposition of two large block circulant matrices. Using fast Fourier transform (FFT) [64], the block circulant tensor decomposition problem in (17) can be re-formulated as following f-diagonal tensor decomposition in the Fourier domain222More details are given in Appendix B.

[TABLE]

Where $\hat{\mathcal{G}}=\mathcal{G}\times_{3}\mathbf{W}$, and $\mathbf{W}$ is the FFT matrix defined as

[TABLE]

The conjugate feature of $\mathbf{W}$ indicates that whichever the threshold technique (low-rank or hard-threshold) is adopted, only the first and second slices in the Fourier domain need to be computed. This strategy can be extended to handle multispectral images by considering only the first $\lfloor N/2\rfloor+1$ slices in the Fourier domain, where $N$ is the number of spectral bands.

III-B Equivalence of f-diagonal Tensor Decomposition in the Fourier Domain

In fact, the computation of (18) does not require explicit formulation of f-diagonal tensor in the Fourier domain. We introduce t-SVD or tensor-SVD [65, 66], a new third-order tensor decomposition framework that demonstrates competitive performance in many applications [67, 68, 69]. t-SVD mainly relies on the definition of t-product $*$ between two third-order tensors using (10) in the Fourier domain.

Definition 1** (t-product).**

The t-product $*$ between two third-order tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and $\mathcal{B}\in\mathbb{R}^{n_{1}\times n_{4}\times n_{3}}$ is an $n_{1}\times n_{4}\times n_{3}$ tensor $\mathcal{C}=\mathcal{A}*\mathcal{B}$ given by

[TABLE]

where $\hat{\mathcal{A}}=\mathcal{A}\times_{3}\mathbf{W}$ and $\hat{\mathcal{B}}=\mathcal{B}\times_{3}\mathbf{W}$.

The t-SVD of a third-order tensor can be defined as the t-product of three third-order tensor via

[TABLE]

which can be computed as

[TABLE]

with an abuse use of t-product, the core tensor can be obtained by

[TABLE]

where $\mathcal{G}*_{i}\mathcal{U}^{T}$ can be computed as the $i$-th mode tensor-matrix product in the Fourier domain. $\mathcal{U}_{group}$ is a third-order tensor with its first frontal slice equal to $\mathbf{U}_{group}$ and other frontal slices equal to $\mathbf{0}$. Furthermore, $\mathcal{U}_{row}$ and $\mathcal{U}_{column}$ can be obtained by minimizing a non-local t-SVD (NL-tSVD) problem

[TABLE]

After some threshold technique, the filtered group $\mathcal{G}_{filtered}$ can be obtained by the inverse transform of (23) via

[TABLE]

III-C Threshold technique

For multiway filtering approaches described in (5), there are roughly two strategies to encourage the sparsity of linear approximation in the transform domain: threshold core tensor $\mathcal{C}$ (via L0-norm [70], Wiener filter [71, 7], soft- or hard-thresholding [7]), and threshold factor matrices via low rank prior [72]. Directly modeling each mode of 4D data with low rank prior raises two major concerns. First, the possible combination of rank estimation along each mode is extensive, and it risks falling into the unbalance trap introduced in [73]. Briefly, given a 4D group of $K$ patches $\mathcal{G}\in\mathbb{R}^{8\times 8\times N\times K}$, its first mode unfolding is a skinny matrix $\mathbf{G}_{1}\in\mathbb{R}^{8\times(8NK)}$, since $8NK\gg 8$, then the rank of $\mathcal{G}$ along the first mode is assumed very low, which risks the loss of more information [73]. To solve these issues, although not explicitly stated, some methods [20, 21] reshape $\mathcal{G}$ into a third-order tensor $\mathcal{G}_{r}\in\mathbb{R}^{64\times N\times K}$ by vectorizing its frontal slices, and impose the low rank constraint only to its grouping dimension. We notice that the block matrix representation in (15) may alleviate the unbalance issue, which may further account for the superiority of low rank t-SVD based methods [69, 74]. However, choosing the multi-rank of f-diagonal tensor (12) in the Fourier domain is not easy. In this paper, we adopt the simple hard-thresholding technique to achieve tensor sparsity.

III-D Global basis and local group representation

The presence of noise in the training process could distort the local representation and introduce some unwanted artifacts. Some recent denoising strategies [75, 76] consider information of the whole image and learn a global representation to render more robustness. CBM3D shares the same idea by applying a pre-defined transform for all patches. t-SVD is also a suitable alternative as a global basis because it preserves the spatial information, and the pre-defined FFT transform along the third mode could make it less sensitive to the variation of noise. The global patch representation can be trained with randomly sampled patches, but for simplicity, all the reference patches are used.

If the patch representation is acquired, then grouping and collaborative filtering could be viewed as a feature extraction and patch classification process that takes care of nonlocal similarity. Therefore, some simple and effective classification method such as PCA can be utilized. We term the combination of global t-SVD basis and local PCA transform as ’multispectral t-SVD (MSt-SVD)’, and detailed implementation is given in Algorithm 1 and Fig. 3. Comparing the FFT matrix $\mathbf{W}$ in (19) and the opponent color mode transform matrix $\mathbf{U}_{color}^{T}$ of CBM3D in (6), it can be seen that the first slice in the Fourier domain corresponding to the first row of $\mathbf{W}$ can be regarded as luminance channel, thus similar to CBM3D, the grouping process and the training of local PCA can be performed by considering only the first slices of all patches in the Fourier domain. Obviously, this implementation can save $\frac{2}{3}$ computational time on grouping and training. This efficient modification of MSt-SVD for color image is termed ’color MSt-SVD’ (CMSt-SVD), and its implementation is briefed in Algorithm 2.

Interestingly, the relationship among nonlocal t-SVD, CBM3D and 4DHOSVD is therefore established using block diagonal representation in the Fourier domain.

III-E Computational complexity

In this subsection, we compare the computational complexity of CBM3D, 4DHOSVD and the proposed MSt-SVD. For simplicity, we assume that the number of image pixels is $N$, that the average time to compute similar patches per reference patch is $T_{s}$, that the average number of patches similar to the reference patch is $K$, and that the size of the patch is $p\times p$ ($p\ll K$). According to [7], the time complexity of 4DHOSVD and CBM3D are $O([T_{s}+Kp^{3}+Kp^{4}]N)$ and $O([T_{s}+Kp^{2}logp+p^{2}KlogK]N)$, respectively. The computational burden of MSt-SVD mainly lies in the PCA transform $O(Kp^{4})$ and patch level t-SVD transform $O(Kp^{3})$, leading to a total complexity of $O([T_{s}+Kp^{3}+Kp^{4}]N)$. Considering that MSt-SVD is a one step algorithm and does not require the training of patch level transform for each group, it is competitive in terms of efficiency.

IV Experiments

In this section, we evaluate the performance of MSt-SVD and CMSt-SVD for color image and multispectral image denoising. All the results of compared methods are obtained by fine-tuned parameters or from the authors’ papers. All the experiments are performed on a moderate laptop equipped with Core(TM) i5-8250U @ 1.8 GHz and 8GB RAM. Our software package is publicly available333https://github.com/ZhaomingKong, which includes a fast C++ mex-function that could help reproduce all our results of color image and multispectral image denoising within 10 and 100 seconds, respectively.

IV-A Experimental setting for color image

A brief description of four publicly available real-world datasets is listed in Table II, and more detailed information is in [43] and [46].

The representative compared methods for color image denoising include: CBM3D [9], 4DHOSVD1 (hard-thresholding) [7], WTR1 [51], Neat Image (NI), TNRD [77], GID [32], MCWNNM [13], TWSC [15], LSCD [78], and LLRT [21]. Three representative neural network based methods MLP [37], DnCNN [38] and FFD-Net [40] are also included in our comparison. Considering the computational complexity of some compared methods and for fair comparison, all methods are tuned to produce their best average results, and the input noise level $\sigma$ of some Gaussian denoisers is listed in Table III. In practical implementation, however, $\sigma$ should be tuned for every image, so to better understand the effectiveness of state-of-the-art CBM3D, the best result of CBM3D on every image is reported, and this implementation is termed ’CBM3D_best’. PSNR and SSIM indices are employed for objective evaluation.

IV-B Experimental results for color image

IV-B1 Experiments on CC15

The fine-tuned results of TNRD, MLP and DnCNN in [13] are used. PSNR result of every image and average computational time are listed in Table IV444CBM3D and LSCD use a C++ mex function, while other methods are implemented purely with Matlab., and visual evaluations are given in Fig. 4 and Fig. 5. Table IV shows that the simple CMSt-SVD consistently outperforms MSt-SVD and 4DHOSVD1, and both MSt-SVD and CMSt-SVD are able to produce very competitive performance in terms of both effectiveness and efficiency. Fig. 4 shows that the representative low-rank based method LLRT and the sparse coding scheme TWSC produce satisfactory results in homogenous regions, because the underlying clean patches share much similar feature, and thus can be modeled as a low-rank or sparse coding problem. But as illustrated in Fig. 5, when the ground truth image contains more details or local variations, clear over-smooth effects can be observed. Interestingly, the state-of-the-art neural network FFD-Net shows similar effects. Besides, compared with CBM3D and CMSt-SVD, the local 4DHOSVD transform is more easily affected by the presence of noise, which is incorporated in the training process of color mode transform.

IV-B2 Experiments on CC60, Xu-100 and our datasets

Table V and VI list the average PSNR and SSIM values of several competitive methods. It is obvious that CMSt-SVD still demonstrates one of the best performance. Our visual evaluations in Fig. 6 and Fig. 8 further illustrate the over-smooth effects of the sparse coding scheme and state-of-the-art neural network, while Fig. 7 shows that many competitive methods including NI, 4DHOSVD, MCWNNM and GID produce color artifacts to some degree, which is similar to our observation in Fig. 4.

IV-B3 Visual evaluation on dnd dataset [42]

The ground truth images of this dataset is not available, so we choose one severely corrupted image that contains lines, smooth regions, color texture and details. The input noise level $\sigma$ of MSt-SVD555An efficient tool is available at github.com/ZhaomingKong/Pure_Image and CBM3D is tuned to achieve the best possible visual effects, and about $\sigma=50$ is used for both MSt-SVD and CBM3D. Also, the parameters of commercial software Neat Image are carefully chosen to compare their difference. Fig. 9 shows the visual comparison. Unfortunately, all compared methods introduce unwanted artifacts. Neat Image presents the best results at line areas, while MSt-SVD produces sharper details with the green color uniformly distributed. The benchmark CBM3D with predefined transforms seeks a balance between details and smoothness.

IV-C Experimental setting for multispectral image

The representative compared methods for multispectral image denoising include: LRTA [19], PARAFAC [79], 4DHOSVD1, LLRT, BM4D [47], TDL [80], ISTReg [81], LRMR [82] and Nmog [35].

Four quality indices are employed for multispectral image: PSNR, SSIM, ERGAS [83] and SAM [84]. EGRAS and SAM are spectral-based evaluation indices, and the smaller EGRAS and SAM values are, the better the restored images are.

Three publicly available datasets are used: Columbia Multispectral Dataset (CAVE)666www1.cs.columbia.edu/CAVE/databases/multispectral for simulated experiments, Harvard Hyperspectral Dataset (HHD)777vision.seas.harvard.edu/hyperspec/index.html [45] and Urban888www.tec.army.mil/hypercube dataset for the real cases.

IV-D Experimental results for Multispectral image

IV-D1 i.i.d. Gaussian Noise

The whole CAVE database consisting of 32 hyperspectral images is used in our synthetic tests. The images of size $512\times 512\times 31$ are captured with the wavelengths in the range of 400-700 nm at an interval of 10 nm. In this experiment, entries in all slices were corrupted by zero-mean i.i.d Gaussian noise $N(0,\sigma^{2})$ with $\sigma=\{10,30,50,100\}$. Since LRMR and ISTReg require much more memory space, their results are copied from [13]. Detailed results are listed in Table VII, visual effect comparison is given in Fig. 10. It can be seen that the recursive use of patch level information by MSt-SVD may better preserve details.

IV-D2 non i.i.d. Gaussian Noise

Entries in all slices were corrupted by zero-mean Gaussian noise with increasing intensity from 21 to 51, and for fair comparison, the input noise level for all compared methods is taken as the average number 36. The detailed results and visual evaluations are given in Table VIII999BM4D uses a C++ mex-function with parallel implementation and TDL uses some mature toolboxes. and Fig. 11, respectively. Comparing Table IV and Table VIII, it is interesting to see that compared with 4DHOSVD, MSt-SVD is faster in dealing with color images but slower in multispectral images, this is because the for-loop of Matlab slice-by-slice computation in the Fourier domain is slow, and our C++ mex-function can reduce the total time to 1.5 minutes.

IV-D3 Experiments on HHD data

There are 50 images of size $1040\times 1392\times 31$, and some of them are clearly contaminated by noise. Considering the large size of noisy images, we mainly examine the effectiveness of MSt-SVD on handling real multispectral image, and compare it with the efficient benchmark TDL. The input noise level $\sigma$ is manually tuned for both TDL and MSt-SVD to balance smooth effects and details, and $N_{step}=8$ is used for MSt-SVD to save some time. Visual evaluation is given in Fig. 12, and the artifacts produced by TDL can be seen in Fig. 12(e) and Fig. 12(k).

IV-D4 Experiments on Urban HSI data

The full Urban dataset of size $307\times 307\times 210$ is used, and some of the bands are contaminated by stripe noise in Fig. 13(a) and Fig. 14(a). As illustrated in 13(b), naively applying MSt-SVD fails to remove stripe noise. In fact, the influence of the sparse stripe noise may be amplified by recursively computing the row- and column-wise relationship. According to (27), the use of block representation (15) does not change the group level representation, thus the influence of stripe noise could be attenuated by adopting sparsity in the grouping dimension. Specifically, we first reshape the original data as a new image of size $210\times 307\times 307$, such that stripe noise sparsely spread along the grouping dimension, and then apply MSt-SVD to the new data, the final filtered image is obtained by reshaping it back to the original size. We term this simple operation ’twist MSt-SVD’ and notice that it does not increase computational burden. Since many compared methods can not handle the stripe noise [35], we use the benchmark Nmog for visual effects comparison. Fig. 13(d) and Fig. 14(c) demonstrate the effectiveness of the twist implementation. Considering its efficiency, its performance may be improved by further modeling of tensor sparsity [70].

IV-E Choice of parameters

Among all free parameters, the hard-threshold parameter $\tau$ directly controls the core tensor sparsity in the transform domain, so we mainly investigate how $\tau$ could influence the proposed MSt-SVD and choose $ps=8$, $SR=20$ and $K$ = 30 according to the settings of 4DHOSVD in [7]. For 4DHOSVD in [7], the authors set $\tau=\sigma\sqrt{2log(n\_elem)}$ based on [85] for their simulated experiments on Kodak gallery101010http://r0k.us/graphics/kodak/, where $n\_elem$ is the number of elements of a patch group. But we observe that it is chosen too large to provide over-smooth effects, so we multiply $\tau$ with a scale factor $\gamma$, and tune the averagely best PSNR value using the first 8 images of Kodak gallery. Fig. 15 shows the influence of $\gamma$ on both 4DHOSVD and MSt-SVD, and the choice of parameters is detailed in Table IX. Fig. 16 gives an illustration of the influence of $\tau$ with Brodatz color texture111111http://multibandtexture.recherche.usherbrooke.ca/colored_brodatz_more.html. In real cases, the tuning of $\tau$ for MSt-SVD is efficient because some intermediate results can be preserved to avoid the recursive computation of grouping and local PCA transform.

V Conclusion

In this paper, we build the relationship among state-of-the-art transforms with block diagonal representation, and investigate the proper choice of patch level transform. According to our discussion and analysis, two simple and effective methods that combine a global t-SVD basis and local PCA transform are proposed. The proposed MSt-SVD and CMSt-SVD utilize more spatial information, and produce competitive performance with state-of-the-art filters in terms of both efficiency and effectiveness. Recently, some statistical properties of tensor decomposition [86] are studied, it is interesting to investigate a further understanding of both color image and multispectral image denoising with block diagonal representation. Besides, our future research also includes classification [87] and related image restoration problems [88].

Appendix A Some Features Related To Block Circulant Representation (15)

Theorem 1**.**

$bcirc(\mathcal{P}_{i})bcirc(\mathcal{P}_{i})^{T}$* is also a block circulant matrix.*

Theorem 2**.**

Two patches $\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ are similar if and only if $bcirc(\mathcal{P}_{i})$ and $bcirc(\mathcal{P}_{j})$ are similar. More specifically,

[TABLE]

Theorem 3**.**

Given a group of similar patches $\mathcal{G}$, and its block circulant tensor representation $bcirc(\mathcal{G})$ in (16). If $\mathbf{U}_{group}$ and $\mathbf{U}_{bcirc\_group}$ are the last mode factor matrix of $\mathcal{G}$ and $bcirc(\mathcal{G})$, respectively, then

[TABLE]

The proof of above Theorem is not hard by checking corresponding definition. Theorem 1 offers an efficient implementation to compute product between a block circulant matrix and its transpose. Theorem 26 and Theorem 27 indicate that nonlocal similarity and linear representation do not change after block circulant operation.

Appendix B Re-formulation of Block Circulant Tensor Decomposition In The Fourier Domain

Theorem 4**.**

[64]** Given a patch $\mathcal{P}_{i}$ and its block circulant representation $bcirc(\mathcal{P})$ in (15), there exists an orthogonal matrix $\mathbf{F}\otimes\mathbf{I}$ that could transform $bcirc(\mathcal{P})$ into a block diagonal matrix via

[TABLE]

where $\otimes$ represents the kronecker product, and $\hat{\mathcal{P}}=\mathcal{P}\times_{3}\mathbf{W}$, $\mathbf{W}=\sqrt{3}\mathbf{F}$. Furthermore, according to the orthogonality, we have

[TABLE]

Acknowledgment

The authors would like to thank all the reviewers for their valuable comments, and appreciate all the authors for sharing their code and software package.