Introduction
噪声水平估计对于非盲去噪方法是至关重要的,噪声水平估计质量直接影响去噪的质量。这篇文章是2015年ICCV 的一篇文章,针对于加性高斯白噪声,其利用非局部相似块具有低秩性的特性,利用协方差矩阵冗余维度的特征值估计噪声水平,并取得了不错的效果。这篇文章的主要贡献在于
- 分析了噪声水平和图像patch的协方差矩阵的特征值之间的关系;
- 提出了一种无参方法在多项式时间内从特征值估计噪声水平的方法,并给与了理论证明;
Analysis
这篇文章的方法是针对于加性,独立,同质性的高斯噪声。所谓“同质性”指的是噪声方差对于图像中所有像素都是一个常数,不会随着位置和颜色强度改变。由于一般认为图像的噪声都是零均值噪声,所谓的噪声水平估计就是通过单张噪声图像估计高斯噪声的方差(或标准差)。
该文章的方法是基于以下观察:从无噪声图像中提取的patch常常处于一个低维子空间中,而不是均匀分布于所有空间中。所以噪声的方差就可以从冗余的空间维度上进行估计,那么噪声估计问题就变成了如何选择冗余维度的问题。
Method
把一张图像分解为一系列的patch
X
s
=
{
x
t
}
t
=
1
s
∈
R
r
×
s
X_s=\{x_t\}^s_{t=1}\in\mathbb{R}^{r\times s}
Xs={xt}t=1s∈Rr×s。如果图像大小为
M
×
N
×
c
M\times N\times c
M×N×c,则
X
s
X_s
Xs包含
s
=
(
M
−
d
+
1
)
(
N
−
d
+
1
)
s=(M-d+1)(N-d+1)
s=(M−d+1)(N−d+1)个大小为
d
×
d
×
c
d\times d\times c
d×d×c的patch。然后将所有这些patch都变成向量,其维度为
r
=
c
d
2
r=cd^2
r=cd2。对于集
X
s
X_s
Xs中任意一个向量
x
t
x_t
xt,其可以分解为:
x
t
=
x
^
t
+
e
t
x_t=\hat{x}_t+e_t
xt=x^t+et
其中, x ^ t \hat{x}_t x^t是对应的无噪声图像patch, e t e_t et表示加性高斯白噪声, e t ∼ N r ( 0 , σ 2 I ) e_t\sim N_r(0,\sigma ^2I) et∼Nr(0,σ2I)。
特征值
从下图中可以看出,无噪声图像patch的特征值多数为0(即说明其位于低维子空间),而噪声图像patch的特征值多数位于真实噪声方差50附近,但不是精确为50。
所以如何精确的获得噪声方差是需要解决的问题。
计算
x
t
x_t
xt的协方差矩阵可得
Σ
x
=
1
s
∑
t
=
1
s
(
x
t
−
μ
)
(
x
t
−
μ
)
T
\Sigma_x=\frac{1}{s}\sum^s_{t=1}(x_t-\mu)(x_t-\mu)^T
Σx=s1t=1∑s(xt−μ)(xt−μ)T
其中,
μ
=
1
s
∑
t
=
1
s
x
t
\mu=\frac{1}{s}\sum^s_{t=1}x_t
μ=s1∑t=1sxt为其平均值。
对其做特征值分解可得
Σ
x
=
R
Λ
R
T
=
R
[
λ
1
0
⋯
0
0
λ
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
λ
r
]
R
T
\Sigma_x=R\Lambda R^T= R\left[\begin{matrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_r \\ \end{matrix}\right]R^T
Σx=RΛRT=R⎣⎢⎢⎢⎡λ10⋮00λ2⋮0⋯⋯⋱⋯00⋮λr⎦⎥⎥⎥⎤RT
则
x
t
x_t
xt可以表示为(此处
y
t
y_t
yt与原文表达意思不同)
x
t
=
R
y
t
x_t=Ry_t
xt=Ryt
假设干净的patch位于
m
m
m维线性空间中,
m
≪
r
m\ll r
m≪r,则
x
^
t
=
A
y
^
t
\hat{x}_t=A\hat{y}_t
x^t=Ay^t
其中,
R
=
[
A
,
U
]
R=[A,U]
R=[A,U],
A
∈
R
r
×
m
A\in \mathbb{R}^{r\times m}
A∈Rr×m表示
m
m
m维的子空间。
综上可知,
x
t
=
R
y
t
=
R
[
y
^
t
,
0
]
+
e
t
=
[
A
,
U
]
[
y
^
t
,
0
]
+
e
t
x_t=Ry_t=R\left[ \begin{matrix} \hat{y}_t,\\ \mathbf{0}\\ \end{matrix}\right]+e_t=\left[ A,U \right] \left[ \begin{matrix} \hat{y}_t,\\ \mathbf{0}\\ \end{matrix}\right]+e_t
xt=Ryt=R[y^t,0]+et=[A,U][y^t,0]+et
公式两边乘以
R
T
R^T
RT
R
T
x
t
=
[
y
^
t
,
0
]
+
R
T
e
t
=
[
y
^
t
+
A
T
e
t
,
U
T
e
t
]
R^Tx_t=\left[ \begin{matrix} \hat{y}_t,\\ \mathbf{0}\\ \end{matrix}\right]+R^Te_t=\left[ \begin{matrix} \hat{y}_t+A^Te_t,\\ U^Te_t\end{matrix}\right]
RTxt=[y^t,0]+RTet=[y^t+ATet,UTet]
通过
R
T
R^T
RT将噪声部分
e
t
e_t
et从
x
t
x_t
xt中分离出来,且根据高斯函数性质,
n
t
=
U
T
e
t
n_t=U^Te_t
nt=UTet是服从高斯分布
N
r
−
m
(
0
,
σ
2
I
)
N_{r-m}(0,\sigma^2I)
Nr−m(0,σ2I)的随机变量。
R
T
x
t
R^Tx_t
RTxt的协方差为
1
s
∑
t
=
1
s
(
R
T
x
t
)
(
R
T
x
t
)
T
=
R
T
Σ
x
R
=
Λ
\frac{1}{s}\sum^s_{t=1}(R^Tx_t)(R^Tx_t)^T=R^T\Sigma_xR=\Lambda
s1t=1∑s(RTxt)(RTxt)T=RTΣxR=Λ
为了估计噪声方差,将特征值
S
=
{
λ
i
}
i
=
1
r
S=\{\lambda_i\}^r_{i=1}
S={λi}i=1r分成两个部分,其中,
S
1
=
{
λ
1
}
i
=
1
m
S_1=\{\lambda_1\}^m_{i=1}
S1={λ1}i=1m表示主维度的特征值,
S
2
=
{
λ
1
}
i
=
m
+
1
r
S_2=\{\lambda_1\}^r_{i=m+1}
S2={λ1}i=m+1r表示冗余维度的特征值。用冗余维度的特征值估计噪声方差。
由于冗余维度向量
n
t
=
U
T
e
t
n_t=U^Te_t
nt=UTet服从高斯分布
N
r
−
m
(
0
,
σ
2
I
)
N_{r-m}(0,\sigma^2I)
Nr−m(0,σ2I),对于向量
n
t
n_t
nt中的每一个元素
n
t
[
i
]
n_t[i]
nt[i]也是服从高斯分布N(0,\sigma^2)的随机变量,由于
λ
i
=
1
s
∑
t
=
1
s
n
t
[
i
]
2
,
i
∈
m
+
1
,
m
+
2
,
.
.
,
r
\lambda_i=\frac{1}{s}\sum^s_{t=1}n_t[i]^2, i\in{m+1,m+2,..,r}
λi=s1t=1∑snt[i]2,i∈m+1,m+2,..,r
则特征值
λ
i
∈
S
2
\lambda_i\in S_2
λi∈S2也是均值为噪声方差
σ
2
\sigma^2
σ2的随机变量。
这里放上文章中的引理:
维度选择
下图表示噪声图像特征值的直方图,可以看出除了少数比较大的异常值(主维度的特征值)以外,其余的特征值是服从高斯分布的,且根据上面的引理,其均值(期望)就是噪声方差。
所以现在的任务就是估计主维度的数量,去除其特征值,即判断特征值是否属于
S
2
S_2
S2。
放上文章另一个定理
其说明,当集合
S
S
S中的冗余维度足够多时,可以通过判断
S
S
S的平均值与其中值是否相等,来判断
S
S
S是否存在异常值。我是从另一个角度理解的,由于冗余维度的特征值服从高斯分布,高斯分布的均值等于其中值,而异常值的存在破坏了这种等价关系,随着异常值的剔除,其中值和均值逐渐接近。
在文章中,patch大小设为
8
×
8
8\times 8
8×8,则
X
=
{
x
t
}
t
=
1
s
X=\{x_t\}^s_{t=1}
X={xt}t=1s维度为192,则主维度的数量不能超过67(具体证明可参见原文)。
算法
根据以上分析,则可以得到噪声水平估计算法,总体来说,就是通过剔除
S
S
S中的异常值,得到冗余维度的特征值,然后冗余维度的特征值的期望(均值)就是噪声的方差。