图像降噪算法——Variance Stabilizing Transform / Generalization Anscombe Transform算法

1. 基本原理

Variance Stabilizing Transform算法指将高斯泊松噪声转化为高斯噪声的一系列算法，对该算法最早的研究开始于1918年，可以说是非常经典了，其在图像降噪上的应用可以参考坦佩雷大学Alessandro Foi组的一系列工作，著名的BM3D算法也是由这个组提出的，这里我们主要介绍其中应用较为广泛的的Generalization Anscombe Transform算法。

首先，我们知道，实际由图像传感器拍摄的图像上的噪声是接近高斯泊松噪声，高斯噪声部分好理解，通常由热噪声造成，而泊松噪声主要是由于光照响应非均匀性造成的，所谓光照响应非均匀性就是说传感器接受光子的不确定性是和信号相关的。下图就展示了高斯泊松噪声的特点，方差曲线随着均值增大而逐渐升高的部分即为泊松噪声，而方差曲线上抖动的部分即为高斯噪声。

*高斯泊松噪声数学建模如下：x^=ap+n\hat{x}=a p+n x^=ap+nη=x^−y^\eta=\hat{x}-\hat{y} η=x^−y^其中η\etaη就是高斯泊松噪声，x^\hat{x}x^为带噪图像，y^\hat{y}y^为无噪图像，aaa为泊松噪声增益，p∼P(y^/a)p \sim \mathcal{P}(\hat{y} /a)p∼P(y^/a)为泊松噪声部分，其均值mmm和方差vvv满足m=v=y^/am=v=\hat{y} /am=v=y^/a**，n∼N(m,σ^2)n \sim \mathcal{N}\left(m, \hat{\sigma}^{2}\right)n∼N(m,σ^2)为高斯噪声部分

我们目前使用的大多数算法，例如NLM, BM3D，都是针对高斯噪声模型设计的，其效果在高斯泊松噪声图像上会大打折扣，如果我们能有这么一种变换将高斯泊松噪声转化为高斯噪声，对高斯噪声降噪后再将图像变换回去，不就非常完美了吗，这种变换就是Generalization Anscombe Transform算法。下面介绍该变换：

正变换：f(x^)={2aax^+38a2+σ^2−am,x^>−38a−σ^2a+m0,x^≤−38a−σ^2a+mf(\hat{x})=\left\{\begin{array}{ll} \frac{2}{a} \sqrt{a \hat{x}+\frac{3}{8} a^{2}+\hat{\sigma}^{2}-a m}, & \hat{x}>-\frac{3}{8} a-\frac{\hat{\sigma}^{2}}{a}+m \\ 0, & \hat{x} \leq-\frac{3}{8} a-\frac{\hat{\sigma}^{2}}{a}+m \end{array}\right. f(x^)={a2ax^+83a2+σ^2−am,0,x^>−83a−aσ^2+mx^≤−83a−aσ^2+m如果我们对增益进行归一化x=x^−ma,σ=σ^ax=\frac{\hat{x}-m}{a}, \sigma=\frac{\hat{\sigma}}{a} x=ax^−m,σ=aσ^也就是说如果我们将输入图像归一化为一个单位泊松变量叠加一个均值为零，标准差为σ\sigmaσ的高斯噪声，那么该变换可以简化为fσ(x)={2x+38+σ2,x>−38−σ20,x≤−38−σ2f_{\sigma}(x)=\left\{\begin{array}{ll} 2 \sqrt{x+\frac{3}{8}+\sigma^{2}}, & x>-\frac{3}{8}-\sigma^{2} \\ 0, & x \leq-\frac{3}{8}-\sigma^{2} \end{array}\right. fσ(x)={2x+83+σ2,0,x>−83−σ2x≤−83−σ2
逆变换：
设yyy为降噪后的图像，闭合形式可以近似为fσ(y)−1=14y2+1432y−1−118y−2+5832y−3−18−σ2f_{\sigma}(y)^{-1}= \frac{1}{4} y^{2}+\frac{1}{4} \sqrt{\frac{3}{2}} y^{-1}-\frac{11}{8} y^{-2}+\frac{5}{8} \sqrt{\frac{3}{2}} y^{-3}-\frac{1}{8}-\sigma^{2} fσ(y)−1=41y2+4123y−1−811y−2+8523y−3−81−σ2当σ\sigmaσ和yyy特别大时，其渐近逆为：fσ(y)−1=14y2−18−σ2f_{\sigma}(y)^{-1}=\frac{1}{4} y^{2}-\frac{1}{8}-\sigma^{2} fσ(y)−1=41y2−81−σ2

我们可以来看下Generalized Anscombe Transform算法正变换后方差如下图(b)所示：

可以看到，当σ=0.01,,1,2,3\sigma=0.01,,1,2,3σ=0.01,,1,2,3时，通过正变换后方差稳定在1左右，而不随yyy发生变化了，也就是通过Generalization Anscombe Transform算法将高斯泊松噪声转变成了高斯噪声

2. matlab代码

这个算法的代码很简单，就不贴在这里了，大家可以去参考matlab代码Denoising software for Poisson and Poisson-Gaussian data，这里有两点需要注意：

Generalization Anscombe Transform如果采用的是归一化形式，在处理图像时需要先对图像进行归一化操作
Generalization Anscombe Transform并没有降噪效果，仅仅是改变了噪声的形态，代码中实际使用的降噪算法还是BM3D或者NLM，具体实验结果大家去看论文就好

3. 补充

Anscombe Transform为Generalization Anscombe Transfrom的特殊形式，Generalization Anscombe Transform是将高斯泊松噪声转换为近似高斯分布，而Anscombe Transform可以将泊松分布转换为近似高斯分布。例如，泊松分布变量x∼P(x^)x \sim \mathcal{P}(\hat{x})x∼P(x^)，有m=v=x^m=v=\hat{x}m=v=x^，Anscombe Transform正变换为：A:x→2x+38A: x \rightarrow 2 \sqrt{x+\frac{3}{8}} A:x→2x+83Anscombe Transform逆变换为：A−1:y→(y2)2−38A^{-1}: y \rightarrow\left(\frac{y}{2}\right)^{2}-\frac{3}{8} A−1:y→(2y)2−83以上形式逆变换为代数逆，该逆变换会给均值引入不友好的偏置，因此还有一个精确无偏逆的闭合形式的近似解：y→14y2+1432y−1−118y−2+5832y−3−18y \rightarrow \frac{1}{4} y^{2}+\frac{1}{4} \sqrt{\frac{3}{2}} y^{-1}-\frac{11}{8} y^{-2}+\frac{5}{8} \sqrt{\frac{3}{2}} y^{-3}-\frac{1}{8} y→41y2+4123y−1−811y−2+8523y−3−81
Genrealization Anscombe Transform的推到过程可以参考博客VST变换
我们可以从另一个角度来理解高斯泊松噪声这个问题，高斯泊松噪声实际上就是信号强的地方噪声大，信号小的地方噪声小这样一个噪声分布，Generalization Anscombe Transform是通过将不同区域噪声变得一致来解决这个问题，那在设计降噪算法的时候是不是也可以通过在不同区域实施不同强度的降噪水平来解决这个问题呢？

此外，这里我写一个各种算法的总结目录图像降噪算法——图像降噪算法总结，对图像降噪算法感兴趣的同学欢迎参考