噪声：Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise

Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise

在文章 CMOS图像传感器中的噪声来源分析中分析了图像传感器中的噪声来源， 噪声：Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data中分析了图像传感器器中的噪声分布近似服从泊松-高斯，即图像传感器中的噪声是一个 异方差的高斯噪声，这一结论对于后续的降噪而言是非常重要的，因为目前应用广泛的降噪算法，诸如：BLF、NLM、BM3D、GuideFilter的一个基础假设是图像当中的噪声为高斯白噪声，很明显泊松-高斯噪声并不符合这一假设，所以出现了对异方差高斯噪声的稳定化变换，称为噪声稳定化(Variance Stabilization Transformation)，本文简称其为VST，即将原本的泊松-高斯噪声转化为高斯白噪声，这一方法由噪声分析及降噪方面的专家 Alessandro Foi提出，并在各ISP中得以广泛应用，参考 Benchmark – Darmstadt Noise Dataset中给出的各降噪算法的客观评价指标，可以看出，经过VST之后的降噪结果的客观评价指标均高于没有经过VST的，足以表明VST对泊松-高斯噪声的降噪影响，关于VST的发展历程可以参考如下几篇论文：

• The transformation of Poisson, binomial and negative binomial data
• A closed-form approximation of the exact unbiased inverse of the Anscombe variance-stabilizing transformation
• On the Inversion of the Anscombe Transformation in Low-Count Poisson Image Denoising
• Optimal inversion of the Anscombe transformation in low-count Poisson image denoising
• Poisson-Gaussian Denoising Using the Exact Unbiased Inverse of the Generalized Anscombe Transformation
• Variance Stabilization for Noisy+Estimate Combination in Iterative Poisson Denoising
• Variance Stabilization in Poisson Image Deblurring

generalized Anscombe transformation

根据泊松-高斯噪声模型，图像的噪声方差可以表示为：
$$V=\alpha y+{\sigma }^2\text{\hspace{1em}(1)}$$
受噪声影响的观测信号可以表示为：
$$z=\alpha p+n\text{\hspace{1em}(2)}$$
其中$z$表示观测信号，$\alpha$为式(1)中的噪声参数，$p$为受泊松噪声影响的信号p?P(y)，E{p}=yp\sim\mathcal{P}\left(y\right)，E\{p\}=yp?P(y)，E{ p}=y，$n$为高斯噪声。式(1)表明图像的噪声强度$V$是随信号$y$而变化的，假设存在某种变换$F$使得经过此变换后图像的噪声方差恒为1，即噪声由泊松-高斯噪声变为高斯白噪声，则可将上式目的公式化：
Var{F(z)}=1(3)Var\{F(z)\}=1\tag{3} Var{ F(z)}=1(3)

• 一维随机变量函数的期望与方差的近似计算方法
  设一维连续随机变量$X$，其期望和方差分别为$\mu 、{\sigma }^2$，令$Y=H(X)$，则随机变量$Y$的期望和方差如下：
  E{Y}≈H(μ)+H′′(μ)2σ2Var{Y}≈H′(μ)2σ2E\{Y\}\approx H(\mu)+\frac{H^{''}(\mu)}{2} \sigma^{2}\\ Var\{Y\}\approx H^{'}(\mu)^{2} \sigma^2 E{ Y}≈H(μ)+2H′′(μ)?σ2Var{ Y}≈H′(μ)2σ2
  以上内容的详细推导可以参考：王新洲.非线性模型参数估计理论与应用

那么式(3)就可以写为：
Var{F(z)}=(dF(z)dz)2Var{z}=1(4)Var\{F(z)\}=(\frac{dF(z)}{dz})^{2}Var\{z\}=1\tag{4} Var{ F(z)}=(dzdF(z)?)2Var{ z}=1(4)
即：
dF(z)dz=1Var{z}(5)\frac{dF(z)}{dz}=\sqrt{\frac{1}{Var\{z\}}}\tag{5} dzdF(z)?=Var{ z}1? ?(5)
对式(5)进行积分就得到了对任意分布噪声的VST，称为generalized Anscombe transformation，将式(1)代入式(5)，并对$z$积分：
$$F(z)=\frac{2}{\alpha }\sqrt{{\sigma }^2+\alpha p}\text{\hspace{1em}(6)}$$
即上式(6)即是对泊松-高斯噪声的VST，同时也可以看出，不仅泊松-高斯分布的噪声能被稳定化为高斯白噪声，只要能显示的求出式(5)的积分，就能完成对相应分布的噪声的稳定化，积分的复杂度直接决定了VST的复杂程度。
上式(6)即是对含泊松-高斯噪声的信号$z$的稳定化过程，将泊松-高斯噪声做如下定义：
$$\eta =z?\alpha y$$
因此噪声分布服从泊松-高斯分布的图像的去噪过程即是从包含泊松-高斯噪声的信号$z$中恢复出原始信号$y$，上述过程中我们假设服从高斯分布的噪声$n$的均值为0，但事实上可能并不是这样的，回顾CMOS图像传感器中的噪声来源分析一文，图像传感器采集下来的Raw数据当中一般包含黑电平，在ISP当中一般直接标定黑电平后从Raw数据中减去黑电平，但这一做法必定会带来噪声信号在接近黑电平附近的点的被减完黑电平之后小于0，随后被Clipping至0，上述Clipping对噪声分布的影响参见文章噪声：Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data，简单来说Clipping过程如下：
$$\tilde{z}=max(z?BLC,0)$$
其中$BLC$为黑电平值，如果想要在降噪前噪声信号不被Clipping的话，我们可以令高斯分布的噪声$n$的均值$u=BLC$，此时可以将式(6)所表示的噪声稳定化过程扩展为：
$$F(z)=\frac{2}{\alpha }\sqrt{{\sigma }^2+\alpha p?\alpha u}\text{\hspace{1em}(7)}$$
最后，参见《Image Processing and Data Analysis》，在上式(7)中引入常数项，扩展式(7)的普适性：
$$F(z)=\frac{2}{\alpha }\sqrt{{\sigma }^2+\alpha p?\alpha u+\frac{3}{8}{\alpha }^2}\text{\hspace{1em}(8)}$$
再由于我们不可能预先获取信号$y$，也不可能将高斯噪声从泊松-高斯噪声中剥离出来，我们只能用含噪信号$z$代替上式中的$p$来近似，则对泊松-高斯噪声的稳定化过程表示为：
$$F(z)=\frac{2}{\alpha }\sqrt{\alpha z+\frac{3}{8}{\alpha }^2+{\sigma }^2?\alpha u}\text{\hspace{1em}(9)}$$
这里值得注意的一点是，VST是一种域变换，其目的是为了稳定噪声，不会影响信号本身的分布形式，即经过VST之后，信号依旧处于空间域当中，只是处于另一个亮度空间当中，在VST域降噪之后需要将降噪结果逆变换为原亮度空间。

Exact Unbiased Inverse Transformation

如果我们已知受泊松-高斯噪声影响的信号$z$，及其相关噪声参数$(\alpha ,\sigma )$，为表示方便，令$z\equiv \frac{z?u}{\alpha }$，则上式噪声稳定化过程简写为：
$$f_{\sigma }(z)=2\sqrt{z+\frac{3}{8}+{\sigma }^2}\text{\hspace{1em}(10)}$$
将噪声稳定化后的含噪信号的降噪结果记为D=E{fσ(z)∣y,σ}D=E\{f_{\sigma}(z)|y,\sigma\}D=E{ fσ?(z)∣y,σ}，降噪后还需要将结果$D$由VST域转换值原本图像域，由于噪声稳定化过程$f_{\sigma }$是非线性的，所以有：
fσ?1(E{fσ(z)∣y})≠E{z∣y}(11)f_{\sigma}^{-1}\left(E\left\{f_{\sigma}(z) \mid y\right\}\right) \neq E\{z \mid y\}\tag{11} fσ?1?(E{ fσ?(z)∣y})??=E{ z∣y}(11)
这意味值直接将信号$D$用$f_{\sigma }^{?1}$进行逆变换将得到原信号$y$的有偏估计，因此需要定义从VST后的信号期望E{fσ(z)∣y}E\left\{f_{\sigma}(z) \mid y\right\}E{ fσ?(z)∣y}到原信号期望E{z∣y}E\{z \mid y\}E{ z∣y}间的映射：
Iσ:E{fσ(z)∣y,σ}?E{z∣y,σ}(12)\mathcal{I}_{\sigma}: E\left\{f_{\sigma}(z) \mid y, \sigma\right\} \longmapsto E\{z \mid y, \sigma\}\tag{12} Iσ?:E{ fσ?(z)∣y,σ}?E{ z∣y,σ}(12)
并且假定E{z∣y,σ}=yE\{z \mid y, \sigma\}=yE{ z∣y,σ}=y，论文《Optimal inversion of the Anscombe transformation in low-count Poisson image denoising》中给出了求解E{fσ(z)∣y,σ}E\left\{f_{\sigma}(z) \mid y, \sigma\right\}E{ fσ?(z)∣y,σ}的过程：
E{fσ(z)∣y,σ}=∫?∞+∞fσ(z)p(z∣y,σ)dz=∫?∞+∞2z+38+σ2∑k=0+∞(yke?yk!2πσ2e?(z?k)22σ2)dz\begin{aligned} E\left\{f_{\sigma}(z) \mid y, \sigma\right\}&=\int_{-\infty}^{+\infty} f_{\sigma}(z) p(z \mid y, \sigma) d z \\ &= \int_{-\infty}^{+\infty} 2 \sqrt{z+\frac{3}{8}+\sigma^{2}} \sum_{k=0}^{+\infty}\left(\frac{y^{k} e^{-y}}{k ! \sqrt{2 \pi \sigma^{2}}} e^{-\frac{(z-k)^{2}}{2 \sigma^{2}}}\right) d z \end{aligned} E{ fσ?(z)∣y,σ}?=∫?∞+∞?fσ?(z)p(z∣y,σ)dz=∫?∞+∞?2z+83?+σ2 ?k=0∑+∞?(k!2πσ2 ?yke?y?e?2σ2(z?k)2?)dz?
下图中c是根据以上积分计算得到的几条逆变换映射曲线。

Low-Count Inverse Transformation

上述对VST的逆变换过程虽然能够精准无偏的重建图像亮度，但计算复杂度过高，我在关于VST的实验过程中，计算一条无偏映射曲线大概需要30分钟，相当缓慢，而且这条曲线本身是$\sigma$的函数，在实际应用过程中也不太可能每个$\sigma$都去计算一条映射曲线，只会是计算几个$\sigma$值所对应的映射曲线，其余$\sigma$值的映射曲线用已知曲线进行插值的方式获取，因此很难做到真正的无偏映射，作者在其余几篇论文中陆续给出了几个近似的VST逆变换公式，总结如下：
$$\begin{array}{rl}{\mathcal{I}}_A(D)& =f^{?1}(D)={\left(\frac{D}{2}\right)}^2?\frac{3}{8}\\ {\mathcal{I}}_B(D)& ={\left(\frac{D}{2}\right)}^2?\frac{1}{8}\\ {\tilde{\mathcal{I}}}_C(D)& =\frac{1}{4}{D}^2+\frac{1}{4}\sqrt{\frac{3}{2}}{D}^{?1}?\frac{11}{8}{D}^{?2}+\frac{5}{8}\sqrt{\frac{3}{2}}{D}^{?3}?\frac{1}{8}\end{array}$$