VST
通常的图像去噪算法常常假设图像的噪声模型为一个加性噪声模型,并且噪声假设为高斯白噪声,即
s
(
x
)
=
s
0
(
x
)
+
n
(
x
)
s(x) = s_0(x)+n(x)
s(x)=s0(x)+n(x)
这里
x
x
x为信号的索引坐标,在图像中常常用二维的坐标来索引
s
0
s_0
s0为原始信号,
s
s
s为观测信号,
n
n
n为噪声,并且
n
∼
N
(
0
,
σ
2
)
n\sim {N(0,\sigma^2)}
n∼N(0,σ2)。如下流程图所示:
然而在现实的物理过程中,有许多需要去噪的过程并不是仅仅被高斯噪声所干扰,针对这样的过程,一类处理方法是对噪声进行重新建模,并且根据新建模的噪声设计特定的去噪过程,另一类处理方法是将新建模的噪声过程转换成高斯白噪声,然后采用已经研究过的针对高斯白噪声有效的去噪算法,两种方法各有优劣,但是第二种方法有两个好处,其一直方便模块化,第二是直接利用现有的去噪算法,避免重新投入资源尽心开发,节约时间和研发成本。这个将特定噪声转换成高斯白噪声的过程我们称之为VST,其反过程是将去噪后的信号转换到原始分布,称之为inverse VST。上述整体流程如下图所示:
针对于CCD/CMOS成像系统,观测信号通常建模为泊松-高斯联合分布,其中高斯噪声对应于感光器件本身的噪声,而泊松过程对应于光子打到感光器件上这样一个计数过程,二者相互独立,构成成像系统整体的噪声建模。如上图中的第一个流程图所示。
泊松-高斯联合过程的建模和参数估计不是本文的重点,稍后会有另一片文档详细介绍其建模过程和参数估计方法,本文的侧重点是在知道了泊松高斯联合分布的参数之后,如何将其变换一个稳定高斯噪声(即VST),以及在去完噪声之后如何通过一个反变换转换到原始信号(即inverse VST).
1.VST过程
泊松-高斯联合分布的建模如下式所示:
s
(
x
)
=
α
p
(
s
0
(
x
)
)
+
n
(
x
)
s(x)=\alpha p(s_0(x))+n(x)
s(x)=αp(s0(x))+n(x)
其中
α
\alpha
α为增益,
n
n
n为均值为
m
m
m,标准差为
σ
\sigma
σ的高斯分布,
p
p
p为依赖于信号的泊松分布,其参数为
λ
0
\lambda_0
λ0.
为了方便,我们在下面的推导过程中省去索引坐标
x
x
x.
我们的目标是寻找一个变换
y
=
f
(
s
)
y=f(s)
y=f(s)使得其方差与原始信号
s
s
s无关。从信号的建模公式可知,其方差为
V
a
r
(
s
)
=
σ
2
+
α
2
λ
0
Var(s)=\sigma^2+\alpha^2\lambda_0
Var(s)=σ2+α2λ0,假设信号
s
s
s变化足够小,在
s
s
s的一个小区域内一次逼近就能够达到足够小的误差,那么
s
s
s与
f
(
s
)
f(s)
f(s)的方差关系为
V
a
r
(
f
)
=
(
d
f
f
s
)
2
V
a
r
(
s
)
Var(f)=(\frac{df}{fs})^2Var(s)
Var(f)=(fsdf)2Var(s),不失一般性,设
V
a
r
(
f
)
=
1
Var(f)=1
Var(f)=1,那么
d
f
d
s
=
1
σ
2
+
α
2
m
0
\frac{df}{ds}=\frac{1}{\sqrt{\sigma^2+\alpha^2m_0}}
dsdf=σ2+α2m01
做一个简单的一阶逼近,
α
m
0
=
x
−
g
\alpha m_0=x-g
αm0=x−g,于是又
d
f
d
s
=
1
α
2
−
α
g
+
α
s
\frac{df}{ds}=\frac{1}{\sqrt{\alpha^2-\alpha g + \alpha s}}
dsdf=α2−αg+αs1
经过简单的变量替换,容易求得上述微分方程的解析解:
f
(
s
)
=
2
α
α
s
+
σ
2
−
α
g
f(s)=\frac{2}{\alpha}\sqrt{\alpha s + \sigma^2-\alpha g}
f(s)=α2αs+σ2−αg
上述推导过程是基于一个简单的线性逼近的假设基础上进行推导的。更一般的,我们希望寻找这样的一个变换,即
y
=
f
(
s
)
=
s
+
c
y=f(s)=\sqrt{s+c}
y=f(s)=s+c.这里需要用到级数展开的一系列理论。具体过程如下:
设
E
(
s
)
=
m
E(s)=m
E(s)=m,令
t
=
s
−
m
t=s-m
t=s−m和
m
′
=
m
+
c
m^{'}=m+c
m′=m+c.对
y
y
y进行级数展开:
y
=
m
′
+
t
=
m
′
[
1
+
1
2
t
m
′
2
−
1
8
t
2
m
′
2
+
1
16
t
3
m
′
3
−
5
128
t
4
m
′
4
+
⋯
]
y=\sqrt{m^{'}+t}=\sqrt{m^{'}}[1+\frac{1}{2} \frac{t}{m^{'2}}-\frac{1}{8} \frac{t^2}{m^{'2}}+\frac{1}{16}\frac{t^3}{m^{'3}}-\frac{5}{128}\frac{t^4}{m{'4}}+\cdots]
y=m′+t=m′[1+21m′2t−81m′2t2+161m′3t3−1285m′4t4+⋯]
因此,
E
(
y
)
=
m
′
+
t
=
m
′
[
1
+
−
1
8
μ
2
m
′
2
+
1
16
μ
3
m
′
3
−
5
128
μ
4
m
′
4
+
⋯
]
E(y)=\sqrt{m^{'}+t}=\sqrt{m^{'}}[1+-\frac{1}{8} \frac{\mu^2}{m^{'2}}+\frac{1}{16}\frac{\mu^3}{m^{'3}}-\frac{5}{128}\frac{\mu^4}{m{'4}}+\cdots]
E(y)=m′+t=m′[1+−81m′2μ2+161m′3μ3−1285m′4μ4+⋯]
这里
μ
i
\mu_i
μi为随机变量
t
t
t的
i
i
i阶中心距
E
2
(
y
)
=
m
′
[
1
−
m
u
2
4
m
′
2
+
μ
3
8
m
′
3
−
4
μ
4
64
m
′
4
+
μ
2
2
64
m
′
4
+
⋯
]
E^2(y)=m^{'}[1-\frac{mu_2}{4m^{'2}}+\frac{\mu_3}{8m^{'3}}-\frac{4\mu_4}{64m^{'4}}+\frac{\mu_2^2}{64m^{'4}}+\cdots]
E2(y)=m′[1−4m′2mu2+8m′3μ3−64m′44μ4+64m′4μ22+⋯]
于是:
V
a
r
(
y
)
=
m
2
4
m
′
−
μ
3
8
m
′
2
−
μ
2
2
−
5
μ
4
64
m
′
3
+
⋯
Var(y)=\frac{m_2}{4m^{'}}-\frac{\mu_3}{8m^{'2}}-\frac{\mu_2^2-5\mu_4}{64m^{'3}}+\cdots
Var(y)=4m′m2−8m′2μ3−64m′3μ22−5μ4+⋯
为了推导上述方差的解析解,有必要对级数中的分子(各阶中心距)和分母进行各自推导
通过对泊松-高斯联合分布的特征函数进行研究,不难推断出
二阶矩
μ
2
=
σ
2
+
α
2
m
0
\mu_2=\sigma^2+\alpha^2m_0
μ2=σ2+α2m0
三阶距
μ
3
=
m
3
−
3
m
1
m
2
+
2
m
1
3
=
α
3
m
0
\mu_3=m_3-3m_1m_2+2m_1^3=\alpha^3m_0
μ3=m3−3m1m2+2m13=α3m0
四阶距
μ
4
=
α
4
m
0
+
3
(
σ
2
+
α
2
m
0
)
2
\mu_4=\alpha^4m_0+3(\sigma^2+\alpha^2m_0)^2
μ4=α4m0+3(σ2+α2m0)2
由于
m
′
=
m
+
c
m^{'}=m+c
m′=m+c,
1
m
′
=
1
m
[
1
−
c
m
+
c
2
m
2
+
⋯
]
\frac{1}{m^{'}}=\frac{1}{m}[1-\frac{c}{m}+\frac{c^2}{m^2}+\cdots]
m′1=m1[1−mc+m2c2+⋯]
1
m
′
2
=
1
m
2
[
1
−
2
c
m
+
3
c
2
m
2
+
⋯
]
\frac{1}{m^{'2}}=\frac{1}{m^2}[1-2\frac{c}{m}+3\frac{c^2}{m^2}+\cdots]
m′21=m21[1−2mc+3m2c2+⋯]
1
m
′
3
=
1
m
3
[
1
−
3
c
m
+
⋯
]
\frac{1}{m^{'3}}=\frac{1}{m^3}[1-3\frac{c}{m}+\cdots]
m′31=m31[1−3mc+⋯]
将上述公式代入方差
V
a
r
(
y
)
Var(y)
Var(y),可以推出
V
=
α
4
+
σ
2
−
α
g
−
c
α
4
m
−
α
2
8
m
+
14
α
2
64
m
+
⋯
V=\frac{\alpha}{4}+\frac{\sigma^2-\alpha g -c\alpha}{4m}-\frac{\alpha^2}{8m}+\frac{14\alpha^2}{64m}+\cdots
V=4α+4mσ2−αg−cα−8mα2+64m14α2+⋯
忽略高阶无穷小量,于是有
V
=
α
4
+
16
(
σ
2
−
α
g
)
−
16
c
α
+
6
σ
2
64
m
V=\frac{\alpha}{4}+\frac{16(\sigma^2-\alpha g)-16c\alpha+6\sigma^2}{64m}
V=4α+64m16(σ2−αg)−16cα+6σ2
为了让
V
V
V与
m
m
m无关,上式中第二项必须为零,于是有
c
=
3
8
α
+
σ
2
−
α
g
α
c=\frac{3}{8}\alpha + \frac{\sigma^2-\alpha g}{\alpha}
c=83α+ασ2−αg
此时方差为
α
/
4
\alpha / 4
α/4,将
c
c
c带入,并归一化得到
t
=
2
α
α
s
+
3
8
α
2
+
σ
2
−
α
g
t=\frac{2}{\alpha}\sqrt{\alpha s + \frac{3}{8}\alpha^2+\sigma^2-\alpha g}
t=α2αs+83α2+σ2−αg
至此,便推出了Poission-Gaussian联合分布的VST变换公式。
参考文献:
- Image Processsing and data analysis
- Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise