论文笔记8：Decomposable Nonlocal Tensor Dictionary Learning for Multispectral Image Denoising

引言

本文结合MSI的两个特性：空间上的非局部相似性和光谱的全局相关性，提出了一种新的用于MSI去噪的张量字典学习模型。一方面，一个典型的自然场景包含了遍布整个空间的相似的局部块(patches)集合，由微观结构的同源聚集组成。通过在这些非局部相似的块之间求平均值，空间噪声有望得到显著缓解。另一方面，MSI包含大量的光谱冗余。也就是说，在不同波段获得的图像总是高度相关的。通过从这些全局相关的光谱信息中提取主成分，期望消除光谱MSI噪声(次要成分)。

在我们的模型中，我们使用一个分组稀疏正则化器来施加相似的MSI块，在它们的稀疏分解中共享相同的字典原子，从而隐式地平均这些块之间的噪声。此外，通过假设空间和光谱上的冗余字典，所提出的张量字典学习模型可以很容易地分解成一系列低秩张量近似问题。每一个问题都对应于一个利用MSIs的光谱相关特性进行的光谱降维模型，并且可以通过一些现成的高阶统计量来解决。因此，可以缓解光谱冗余问题。

N阶张量	$\mathcal{A}\in {\mathbb{R}}^{I_1\times \cdots \times I_n\times \cdots I_N}$
$\mathcal{A}$中元素	$a_{i_1\cdots i_n\cdots i_N},1⩽i_n⩽I_n$
n模向量	对$a$改变指标$i_n$固定其他指标
n模向量作为列构成的矩阵	${\mathbf{A}}_{(n)}\in {\mathbb{R}}^{I_n\times (I_1\cdots I_{n-1}{I}_{n+1}\cdots I_N)}$
$\mathcal{A}$的$n$秩	$r_n$，$\mathcal{A}$的n模向量span的向量空间的维数
n模乘积	$\mathcal{A}{\times }_n\mathbf{B}(\in {\mathbb{R}}^{J_n\times I_n})=\mathcal{C}\in {\mathbb{R}}^{I_1\times \cdots \times J_n\times \cdots I_N}$或${\mathbf{C}}_{(n)}=\mathbf{B}{\mathbf{A}}_{(n)}$ 见下图
	$c_{i_1\cdots i_{n-1}{j}_n{i}_{n+1}\cdots i_N}={\sum }_{i_n}{a}_{i_1\cdots i_{n-1}{i}_n{i}_{n+1}\cdots i_N}{b}_{j_n{i}_n}$
张量的F范数	见下式

$$\Vert \mathcal{A}{\Vert }_F={\left(\sum _{i_1,\dots ,i_N}{\left|a_{i_1\cdots i_N}\right|}^2\right)}^{1/2}$$

相关工作

• 2D extended approach：NLM、K-SVD、BM3D、BM4D
• Tensor-based approach：张量分解作为多向滤波的一个特例，可以看作是传统奇异值分解(SVD)的推广。LRTA、并行因子分析法等，然而，它们没有利用MSI的非局部相似性，其性能可能对噪声的程度和类型敏感。

LRTA：Denoising and dimensionality reduction using multilinear tools for hyperspectral images
PARAFAC：Denoising of hyperspectral images using the parafac model and statistical performance analysis

可分解的非局部MSI字典学习

From Image DL to MSI DL

对于一组图像块(按字典顺序排列为列向量) { x i } i = 1 n ⊂ R d \left\{\mathbf{x}_{i}\right\}_{i=1}^{n} \subset \mathbb{R}^{d} {xi​}i=1n​⊂Rd，其中$d$为维数，$n$为图像块的个数，字典学习的目标是要计算字典$\mathbf{D}=[{\mathbf{d}}_1,\cdots ,{\mathbf{d}}_m]\in {\mathbb{R}}^{d\times m}$，其中${\mathbf{d}}_i(m>d)$为atoms，系数矩阵为$\mathbf{Z}=[{\mathbf{z}}_1,\cdots ,{\mathbf{z}}_n]\in {\mathbb{R}}^{m\times n}$，由${\mathbf{x}}_i$的表示系数${\mathbf{z}}_i$组成，通过解决以下优化模型得到：$$\underset{\mathbf{D},{\mathbf{z}}_1,\cdots ,{\mathbf{z}}_n}{\min }\sum _{i=1}^n\Vert {\mathbf{x}}_i-\mathbf{D}{\mathbf{z}}_i\Vert \text{ s.t. }\mathcal{P}\left({\mathbf{z}}_i\right)\le k(1)$$其中$\mathcal{P}(\cdot )$为稀疏控制算子，如$l_0$和$l_1$范数。

推广到MSI中来。MSI：三阶张量$\mathcal{H}\in {\mathbb{R}}^{d_W\times d_H\times d_S}$，通过扫描所有重叠的MSI，我们可以从MSI构建一组3D全波段patches(FBP)： { P i , j } 1 ⩽ i ⩽ d W − d w + 1 , 1 ⩽ j ⩽ d H − d h + 1 ⊂ R d w × d h × d S ( d w < d W , d h < d H ) \{\mathcal{P}_{i,j}\}_{1\leqslant i\leqslant d_W-d_w+1,1\leqslant j\leqslant d_H-d_h+1}\subset\mathbb{R}^{d_w\times d_h\times d_S}(d_w<d_W,d_h<d_H) {Pi,j​}1⩽i⩽dW​−dw​+1,1⩽j⩽dH​−dh​+1​⊂Rdw​×dh​×dS​(dw​<dW​,dh​<dH​)，我们将所有的FBP重构为一组三维patches { X i } i = 1 n \{\mathcal{X}_i\}_{i=1}^n {Xi​}i=1n​，其中$n=(d_W-d_w+1)(d_H-d_h+1)$为patch个数。如此构造的每个FBP包含局部空间和全局光谱维度，这可以很容易地帮助我们考虑MSI背后的两个重要属性：空间patches之间的非局部相似性和所有波段之间的全局相关性。

构建MSI字典模型来计算空间和光谱字典 { D W ∈ R d w × m W , D H ∈ R d h × m H , D S ∈ R d S × m S } \{\mathbf{D}^W\in\mathbb{R}^{d_w\times m_W},\mathbf{D}^H\in\mathbb{R}^{d_h\times m_H},\mathbf{D}^S\in\mathbb{R}^{d_S\times m_S}\} {DW∈Rdw​×mW​,DH∈Rdh​×mH​,DS∈RdS​×mS​}，其中$m_W>d_w,m_H>d_h,m_S>d_S$，表明字典的冗余性
$$\begin{array}{c}{\min }_{{\mathbf{D}}^W,{\mathbf{D}}^H,{\mathbf{D}}^S,{\mathcal{Z}}_i}{\sum }_{i=1}^n\Vert {\mathcal{X}}_i-{\mathcal{Z}}_i{\times }_1{\mathbf{D}}^W{\times }_2{\mathbf{D}}^H{\times }_3{\mathbf{D}}^S\mid \\ \text{s.t.},\mathcal{P}\left({\mathcal{Z}}_i\right)\le k(2)\end{array}$$

其中${\mathcal{Z}}_i\in {\mathbb{R}}^{m_W\times m_H\times m_S}$对应于${\mathcal{X}}_i$的系数张量，管理字典之间的关联交互，$\mathcal{P}(\cdot )$为稀疏正则化项如$l_0$或$l_1$算子。

From Image Group-Sparsity to MSI Group-Block-Sparsity

将相似的patches聚成一类， X k = { x i j k } j = 1 n k , k = 1 , 2 , ⋯   , K \mathbf{X}^{k}=\{\mathbf{x}_{i_j^k}\}_{j=1}^{n_k},k=1,2,\cdots,K Xk={xijk​​}j=1nk​​,k=1,2,⋯,K，$K$为类的个数，$n_k$表示第$k$类中的patch的个数，$i_j^k$表示在第$k$个类第$j$个块中的指标，这使得每个类在字典中共享相似的atoms. 第$k$个类${\mathbf{X}}^k$对应的系数矩阵为${\mathbf{Z}}^k=[{\mathbf{z}}_{i_1^k},{\mathbf{z}}_{i_2^k},\cdots ,{\mathbf{z}}_{i_{n_k}^k}]\in {\mathbb{R}}^{m\times n_k}$，这个同时稀疏编码目标可以通过在每个${\mathbf{Z}}^{(k)}$上应用(1)和下面的组稀疏正则化来实现：$${\Vert {\mathbf{Z}}^{(k)}\Vert }_{p,q}=\sum _{i=1}^m{\Vert {\hat{\mathbf{z}}}_i^k\Vert }_q^p$$

其中${\mathbf{z}}_i^k$为${\mathbf{Z}}^{(k)}$的第$i$行，$(p,q)$通常为$(1,2)$或$(0,\infty )$. 这种组稀疏性正则化器有助于施加${\mathbf{Z}}^{(k)}$的一些全零行。

这种非局部方法可以很容易地推广到MSI中，首先，将相似的FBPs聚成一类 { X i j k } j = 1 n k ( k = 1 , 2 , ⋯   , K ) \{\mathcal{X}_{i_j^k}\}_{j=1}^{n_k}(k=1,2,\cdots,K) {Xijk​​}j=1nk​​(k=1,2,⋯,K)，其中$K$为类的个数，$n_k$为第$k$类中FBP的个数，$i_j^k$为第$k$类第$j$块的指标，我们试图强制每个类在每个空间字典$D^W,D^H$和光谱字典$D^S$中共享相似的atoms. 为了方便，我们联合第$k$类的FBPs，形成4阶张量${\mathcal{X}}^{(k)}\in {\mathbb{R}}^{d_w\times d_h\times d_S\times n_k}$，同样第$k$个FBP类的系数张量 { Z i j k } j = 1 n k \{\mathcal{Z}_{i_j^k}\}_{j=1}^{n_k} {Zijk​​}j=1nk​​形成${\mathcal{Z}}^{(k)}\in {\mathbb{R}}^{m_W\times m_H\times m_S\times n_k}$. 然后，通过下面的Group-block-sparsity正则化可以达到非局部MSI张量DL的目的。

看图2就能很容易理解上面的定义，该概念已经被证实能够增强原始高阶信号的更好恢复，因为它隐含地结合了关于真实信号的有价值的先验信息，并且便于在信号表示中充分利用每个模式的字典原子。

那么我们可以构造如下非局部MSI DL模型：
$$\begin{array}{c}{\min }_{{\mathbf{D}}^W,{\mathbf{D}}^H,{\mathbf{D}}^S,{\mathbf{z}}^{(k)}}{\sum }_{k=1}^K\Vert {\mathcal{X}}^{(k)}-{\mathcal{Z}}^{(k)}{\times }_1{\mathbf{D}}^W{\times }_2{\mathbf{D}}^H{\times }_3{\mathbf{D}}^S\Vert \\ \text{s.t.},{\Vert {\mathcal{Z}}^{(k)}\Vert }_B⪯\left(r_k^W,r_k^H,r_k^S\right)\end{array}(4)$$

其中${\mathbf{v}}_1⪯{\mathbf{v}}_2$表示${\mathbf{v}}_1$每个元素小于等于${\mathbf{v}}_2$每个元素，${\mathcal{Z}}^{(k)}$的Group-block-sparsity保证了每个类${\mathcal{X}}^{(k)}$共享字典$D^W,D^H,D^S$的$r_k^W,r_k^H,r_k^S$个atoms.

问题：如何生成FBPs的类以及如何设置Group-block-sparsity阈值$r_k^W,r_k^H,r_k^S$，对于第一个问题，使用k-means++算法获得所有FBPs的类，第二个问题下面来讨论。

Decomposable Nonlocal MSI DL Model

假设字典$D^W,D^H,D^S$足够冗余，使得不同类中使用的字典原子没有重叠。也就是说，我们假设空间和光谱字典可以写成$D^W=[D_1^W,\cdots ,D_K^W],D^H=[D_1^H,\cdots ,D_K^H],D^S=[D_1^S,\cdots ,D_K^S]$，其中$D_k^W\in {\mathbb{R}}^{d_w\times r_k^W},D_k^H\in {\mathbb{R}}^{d_h\times r_k^H},D_k^S\in {\mathbb{R}}^{d_S\times r_k^S}$且${\sum }_{k=1}^K{r}_k^W=m_W,{\sum }_{k=1}^K{r}_k^H=m_H,{\sum }_{k=1}^K{r}_k^S=m_S$，使得每一个MSI类${\mathcal{X}}^{(k)}$只与子字典$D_k^W,D_k^H,D_k^S$有关，这一假设的合理性在于空间和光谱字典的冗余设置，即使当我们假设两个类共享字典的一个原子时，这一假设仍然成立，因为很容易在字典中复制这个原子。在这种假设下，等式(4)可以等效地重新表述为：
$$\begin{array}{l}\Vert {\mathcal{X}}^{(k)}-{\mathcal{Z}}^{(k)}{\times }_1{\mathbf{D}}^W{\times }_2{\mathbf{D}}^H{\times }_3{\mathbf{D}}^S\Vert \\ =\Vert {\mathcal{X}}^{(k)}-Sub\left({\mathcal{Z}}^{(k)}\right){\times }_1{\mathbf{D}}_k^W{\times }_2{\mathbf{D}}_k^H{\times }_3{\mathbf{D}}_k^S\Vert \end{array}(5)$$

其中$Sub\left({\mathcal{Z}}^{(k)}\right)\in r_k^W\times r_k^H\times r_k^S\times n_k$，原来的非局部MSI DL问题可以分解为一系列子问题：
$$\underset{{\mathbf{D}}_k^W,{\mathbf{D}}_k^H,{\mathbf{D}}_k^S,\mathcal{Y}}{\min }\Vert {\mathcal{X}}^{(k)}-\mathcal{Y}{\times }_1{\mathbf{D}}_k^W{\times }_2{\mathbf{D}}_k^H{\times }_3{\mathbf{D}}_k^S\Vert (6)$$

经过这样的转换后，原来有约束的问题(4)现在被重新表述为一系列没有约束的更小的问题。应该注意的是，每个MSI类张量${\mathcal{X}}^{(k)}$在其第三光谱模式中具有维数冗余，这是由于其重要的内在属性之一：跨谱的全局相关性(global correlation across spectrum)。这意味着${\mathcal{X}}^{(k)}$可以由低秩张量近似获得：
$$\underset{{\mathbf{U}}_1,{\mathbf{U}}_2,{\mathbf{U}}_3,{\mathbf{U}}_4,\mathcal{G}}{\min }\Vert {\mathcal{X}}^{(k)}-\mathcal{G}{\times }_1{\mathbf{U}}_1{\times }_2{\mathbf{U}}_2{\times }_3{\mathbf{U}}_3{\times }_4{\mathbf{U}}_4\Vert (7)$$

其中${\mathbf{U}}_1\in {\mathbb{R}}^{d_k^W\times r_k^W},{\mathbf{U}}_2\in {\mathbb{R}}^{d_k^H\times r_k^H},{\mathbf{U}}_3\in {\mathbb{R}}^{d_k^S\times r_k^S},{\mathbf{U}}_4\in {\mathbb{R}}^{d_k^N\times r_k^N}$对应于${\mathcal{X}}^{(k)}$的四个模式的基向量，且$d_k^W\ge r_k^W,d_k^H\ge r_k^H,d_k^S\ge r_k^S,d_k^N\ge r_k^N$. 这里$\mathcal{G}\in {\mathbb{R}}^{r_k^W\times r_k^H\times r_k^S\times r_k^N}$是core tensor，$r_k^S<d_k^S$导致${\mathcal{X}}^{(k)}$的光谱模的维数降低。(7)很容易利用Tucker分解技术求得，然后令${\mathbf{D}}_k^W={\mathbf{U}}_1,{\mathbf{D}}_k^H={\mathbf{U}}_2,{\mathbf{D}}_k^S={\mathbf{U}}_3,\mathcal{Y}=\mathcal{G}{\times }_4{\mathbf{U}}_4$得到(6)的解。

L. R. Tucker. Some mathematical notes on three-mode factor analysis.

对于(7)中秩参数$r_k^W,r_k^H,r_k^S,r_k^N$的选取，我们可以利用著名的AIC/MDL方法对每个类张量${\mathcal{X}}^{(k)}$的模$i(i=1,2,3,4)$展平${\mathbf{X}}_{(i)}^{(k)}$.

M. Wax and T. Kailath. Detection of signals by information theoretic criteria.

Decomposable Nonlocal MSI DL Algorithm

利用${\mathcal{Z}}^{(k)},{\mathbf{D}}^W,{\mathbf{D}}^H,{\mathbf{D}}^S$来恢复所有重叠的FBP，并对结果进行平均以获得去噪的MSI。第二步使用的k-means++，第三步使用的AIC/DIC和第四步使用的Tucker因子分解均可以快速执行，保证了算法的高效性。

实验

补充材料

数据集：http://www1.cs.columbia.edu/CAVE/databases/multispectral

噪声模型：AWGN和Poisson noise，第一组实验我们用不同$sigma$(高达0.3)的高斯噪声和固定$\kappa =5$的泊松噪声来扰动32个哥伦比亚MSI中的每一个。在第二种情况下，我们使用$\kappa$从2到6，固定$\sigma =0.1$.

在实施去噪方法之前，首先通过方差稳定变换(VST)来重构有噪声的MSI，并且在去噪之后，使用相应的逆变换来获得最终的MSI重构。

FBP宽度$d_w$和高度$d_h$是我们的算法中唯一需要设置的两个参数，在我们所有的实验中，我们只是简单地将它们设置为$d_w=d_h=8$.

对比方法：K-SVD，BM3D，3D-cube K-SVD，ANLM3D，BM4D，LRTA，PARAFAC.

评估指标：峰值信噪比(PSNR)、结构相似性(SSIM)、特征相似性(FSIM)、相对全局平均误差(ERGRAS)和光谱角度映射器(SAM)。



自然场景：https://personalpages.manchester.ac.uk/staff/david.foster/Hyperspectral_images_of_natural_scenes_02

可以很容易地观察到，通过我们的方法恢复的图像适当地去除了噪声，同时精细地保留了图像下面的结构，而通过大多数其他竞争方法获得的结果与原始图像相比包含明显的模糊区域。在这些方法中，ANLM3D和LRTA在结构保持方面表现较好。然而，由它们恢复的图像仍然比由我们的方法获得的图像有更多意想不到的尖锐噪声。

源代码

地址：http://gr.xjtu.edu.cn/c/document_library/get_file?folderId=1766524&name=DLFE-38410.zip
除了有张量字典学习的代码，还有KSVD、NLM、BM3D、BM4D、LRTA、PARAFAC等代码，以及PSNR、SSIM、FSIM、ERGAS、SAM代码。

Usage:

For MSI with Gaussian noise, you can simply use the function TensorDL as follows:
    [ clean_img, basis, idx, cores, info ] = TensorDL( noisy_img, params )
Please type 'help TensorDL' to get more information.

For MSI with Poisson noise or mixed Gaussian and Poisson noise, you can follow these steps:
  1. Apply VST via function GenAnscombe_forward with estimated Gaussian STD \sigma.
     After VST, the STD of noise \sigma_VST is 1.
  2. Re-arrange the MSI into [0, 1]. Remember to scale \sigma_VST as well. 
  3. Apply TensorDL as in the Gaussian noise case.
  4. Reverse the operation in step 2 to re-arrange the filtered MSI.
  5. Apply inverse VST by function GenAnscombe_inverse_exact_unbiased (or alternatives provided
     by toolbox invsansc) to get the final result.
You may find example codes in file tdl_demo\Comparison.m

Also, you can use the demo to see some comparison. You can:
  1. Type 'TDLdemo' to see the pre-computed results.
  2. Type 'Comparison' to run various methods (NOTE: this may take a while, since some of the selected
     methods are time-consuming).
  3. Use 'help Comparison' for more information.
  4. Change test MSI by simply modifying variable 'filename' in Comparison.m (NOTE: make sure your MSI
     meets the format requirements).
  5. Change noise level by modifying variables 'kappa' and 'sigma_ratio' in Comparison.m
  6. Select competing methods by turn on/off the enable-bits in Comparison.m