VST
通常的图像去噪算法常常假设图像的噪声模型为一个加性噪声模型，并且噪声假设为高斯白噪声，即
s
(
x
)
=
s
0
(
x
)
+
n
(
x
)
s(x) = s_0(x)+n(x)
s
(
x
)
=
s
0
(
x
)
+
n
(
x
)
这里
x
x
x
为信号的索引坐标，在图像中常常用二维的坐标来索引
s
0
s_0
s
0
为原始信号，
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
(
s
0
(
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
V
a
r
(
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)
V
a
r
(
f
)
=
(
f
s
d
f
)
2
V
a
r
(
s
)
，不失一般性，设
V
a
r
(
f
)
=
1
Var(f)=1
V
a
r
(
f
)
=
1
，那么
d
f
d
s
=
1
σ
2
+
α
2
m
0
\frac{df}{ds}=\frac{1}{\sqrt{\sigma^2+\alpha^2m_0}}
d
s
d
f
=
σ
2
+
α
2
m
0
1
做一个简单的一阶逼近，
α
m
0
=
x
−
g
\alpha m_0=x-g
α
m
0
=
x
−
g
，于是又
d
f
d
s
=
1
α
2
−
α
g
+
α
s
\frac{df}{ds}=\frac{1}{\sqrt{\alpha^2-\alpha g + \alpha s}}
d
s
d
f
=
α
2
−
α
g
+
α
s
1
经过简单的变量替换，容易求得上述微分方程的解析解：
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
+
2
1
m
′
2
t
−
8
1
m
′
2
t
2
+
1
6
1
m
′
3
t
3
−
1
2
8
5
m
′
4
t
4
+
⋯
]
因此，
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
+
−
8
1
m
′
2
μ
2
+
1
6
1
m
′
3
μ
3
−
1
2
8
5
m
′
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]
E
2
(
y
)
=
m
′
[
1
−
4
m
′
2
m
u
2
+
8
m
′
3
μ
3
−
6
4
m
′
4
4
μ
4
+
6
4
m
′
4
μ
2
2
+
⋯
]
于是：
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
V
a
r
(
y
)
=
4
m
′
m
2
−
8
m
′
2
μ
3
−
6
4
m
′
3
μ
2
2
−
5
μ
4
+
⋯
为了推导上述方差的解析解，有必要对级数中的分子（各阶中心距）和分母进行各自推导
通过对泊松-高斯联合分布的特征函数进行研究，不难推断出
二阶矩
μ
2
=
σ
2
+
α
2
m
0
\mu_2=\sigma^2+\alpha^2m_0
μ
2
=
σ
2
+
α
2
m
0
三阶距
μ
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
=
m
3
−
3
m
1
m
2
+
2
m
1
3
=
α
3
m
0
四阶距
μ
4
=
α
4
m
0
+
3
(
σ
2
+
α
2
m
0
)
2
\mu_4=\alpha^4m_0+3(\sigma^2+\alpha^2m_0)^2
μ
4
=
α
4
m
0
+
3
(
σ
2
+
α
2
m
0
)
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
=
m
1
[
1
−
m
c
+
m
2
c
2
+
⋯
]
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
′
2
1
=
m
2
1
[
1
−
2
m
c
+
3
m
2
c
2
+
⋯
]
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
′
3
1
=
m
3
1
[
1
−
3
m
c
+
⋯
]
将上述公式代入方差
V
a
r
(
y
)
Var(y)
V
a
r
(
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
α
+
4
m
σ
2
−
α
g
−
c
α
−
8
m
α
2
+
6
4
m
1
4
α
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
α
+
6
4
m
1
6
(
σ
2
−
α
g
)
−
1
6
c
α
+
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
=
8
3
α
+
α
σ
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
+
8
3
α
2
+
σ
2
−
α
g
至此，便推出了Poission-Gaussian联合分布的VST变换公式。
参考文献：
- Image Processsing and data analysis
- Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise
