噪声：Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data

Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data

Sensor的噪声特性是确定一款特定Sensor是否满足最终产品规格需求的重要指标，这篇文章是Alessandro Foi于2007年发表的论文《 Noise Measurement for Raw-Data of Digital Imaging Sensors by Automatic Segmentation of Nonuniform Targets》的基础上所提出的用于Sensor噪声分析的框架，其噪声模型也是目前市场上常见的芯片噪声分析(Noise Proflie)模块的基础模型，Alessandro Foi是芬兰Tampere大学研究图像去噪和噪声建模的教授，他的文章谈到的方法很多都已经在手机camera中得到了应用。不同于Imatest所使用到的基于傅里叶变换的方法 ISO 15739 — Noise measurements，本文所提的方法基于简单的图像分割，简单易行。
随着Sensor制成工艺的日益进步，Sensor的像素尺寸越来越小，其所能够积累的光子数量越来越少，最终造成Sensor采集下来信号(Raw)信噪比的下降，同时Sensor的噪声也并不是一般图像处理过程中所假设的高斯白噪声，因此对Sensor所采集图像信号的噪声的正确建模对于后续的图像处理过程是十分重要的。

Poissonian-Gaussian Modeling

参考《CMOS图像传感器中的噪声来源分析》一文中对CMOS Sensor的噪声来源的分析，我们知道Sensor采集下来的Raw数据当中最主要的噪声是泊松噪声和高斯噪声，且这两部分噪声都是白噪声，将这两部分噪声的观测模型做如下表示：
$$z(x)=y(x)+\sigma (y(x))\xi (x)\text{\hspace{1em}(1)}$$
其中$x\in X$表示像素坐标，$z$为当前点被噪声影响后的观测像素值，$y$为当前点不被噪声影响的实际像素值，$\xi (x)$表示均值为0，标准差为1的随机变量，$\sigma$为$y$的函数，表示噪声的标准差随实际像素值变化，$\sigma (y(x))\xi (x)$表示Raw数据的噪声是均值为0、方差随实际像素值变化的而变化随机变量，即可以将Sensor的噪声模型看做是一个异方差的高斯噪声( Heteroskedastic Normal)。在接下来的描述当中将使用 E ⁡ { ⋅ } \operatorname{E}\{\cdot\} E{ ⋅}表示随机变量的数学期望， var ⁡ { ⋅ } 、 std ⁡ { ⋅ } \operatorname{var}\{\cdot\}、\operatorname{std}\{\cdot\} var{ ⋅}、std{ ⋅}分别表示随机变量的方差与标准差。由于 E { ξ ( x ) } = 0 E\{\xi(x)\}=0 E{ ξ(x)}=0，就有 E { z ( x ) } = y ( x ) E\{z(x)\}=y(x) E{ z(x)}=y(x)，原始信号可以被看做是噪声信号的期望，而噪声信号的方差与观察信号、原始信号间的关系如下：
s t d { z ( x ) } = σ { y ( x ) } = σ ( E { z ( x ) } ) (2) std\{z(x)\}=\sigma\{y(x)\}=\sigma(E\{z(x)\})\tag{2} std{ z(x)}=σ{ y(x)}=σ(E{ z(x)})(2)
Sensor当中的泊松噪声来源于其散粒噪声及暗电流，高斯噪声来源于热噪声、复位噪声、量化噪声，可以将这两类噪声看做是相互独立的两部分，分别用${\eta }_p、{\eta }_g$表示如下：
$$\sigma (y(x))\xi (x)={\eta }_{\mathrm{p}}(y(x))+{\eta }_{\mathrm{g}}(x)\text{\hspace{1em}(3)}$$

$$\chi \left(y(x)+{\eta }_{\mathrm{p}}(y(x))\right)\sim \mathcal{P}(\chi y(x)),{\eta }_{\mathrm{g}}(x)\sim \mathcal{N}(0,b)\text{\hspace{1em}(4)}$$

其中$\mathcal{P}$表示泊松分布，$\mathcal{N}$表示高斯分布，$\chi ＞0、b\ge 0$分别为尺度参数。由于泊松分布的方差与均值相等，即有：
E { χ ( y ( x ) + η p ( y ( x ) ) ) } = var ⁡ { χ ( y ( x ) + η p ( y ( x ) ) ) } = χ y ( x ) (5) E\left\{\chi\left(y(x)+\eta_{\mathrm{p}}(y(x))\right)\right\}=\operatorname{var}\left\{\chi\left(y(x)+\eta_{\mathrm{p}}(y(x))\right)\right\}=\chi y(x)\tag{5} E{ χ(y(x)+ηp​(y(x)))}=var{ χ(y(x)+ηp​(y(x)))}=χy(x)(5)
又由于 var ⁡ { χ ( y ( x ) + η p ( y ( x ) ) ) } = χ 2 v a r { η p ( y ( x ) ) } = χ y ( x ) \operatorname{var}\left\{\chi\left(y(x)+\eta_{\mathrm{p}}(y(x))\right)\right\}=\chi ^{2}var\{\eta_{p}(y(x))\}=\chi y(x) var{ χ(y(x)+ηp​(y(x)))}=χ2var{ ηp​(y(x))}=χy(x)，即有 var ⁡ { η p ( y ( x ) ) } = y ( x ) / χ \operatorname{var}\left\{\eta_{\mathrm{p}}(y(x))\right\}=y(x) / \chi var{ ηp​(y(x))}=y(x)/χ且 E { η p ( y ( x ) ) } = 0 E\left\{\eta_{\mathrm{p}}(y(x))\right\}=0 E{ ηp​(y(x))}=0，因此泊松噪声的方差随$y(x)$而变化，记为：
var ⁡ { η p ( y ( x ) ) } = a y ( x ) ; a = χ − 1 (6) \operatorname{var}\left\{\eta_{\mathrm{p}}(y(x))\right\}=a y(x);a=\chi^{-1}\tag{6} var{ ηp​(y(x))}=ay(x);a=χ−1(6)
因此式(1)中的噪声部分可以做如下表示：
$${\sigma }^2(y(x))=ay(x)+b\text{\hspace{1em}(7)}$$

作：原图；右：添加a=0.01、b=0.0016泊松-高斯噪声后的观测图

上图为向原图$y$中所添加噪声的噪声参数$\alpha 、b$取不同值时，最后通过本文方法标定出来的信号期望与噪声方差间的关系示意图，其中实线为添加噪声时所假设的噪声模型，对应的虚线是本文方法对该噪声的标定结果，在首尾两端出现的不一致是由于收尾两端的像素值受噪声干扰后很容易出现小于0或者大于255的情况，最后分别被Clip至0或255造成的(在图像处理中一般需要做Clip的原因是，希望被处理图像数据的灰度值范围固定在8bit、10bit、12bit的范围内，以节省计算、存储等资源)。
为了更进一步的简化噪声模型，用高斯分布来进一步近似泊松分布，即：
$$\mathcal{P}(\lambda )\approx \mathcal{N}(\lambda ,\lambda )$$
当$\lambda$足够大时，可以将泊松分布看做是异方差的高斯分布，此时泊松-高斯分布的噪声模型可以进一步简化为：
$$\sigma (y(x))\xi (x)=\sqrt{ay(x)+b}\xi (x)\simeq {\eta }_{\mathrm{h}}(y(x))\text{\hspace{1em}(8)}$$
其中${\eta }_{\mathrm{h}}(x)\sim \mathcal{N}(0,ay(x)+b)$。

The Noise Profile Algorithm

Wavelet domain analysis

通常情况下，我们并不能直接获取$(\sigma ,y)$的数据对来拟合式(8)所示噪声模型的参数$\alpha ,b$，用小波分解的方式将图像分解为低频部分$z^{wapp}$、高频部分$z^{wdet}$，其中低频部分包含了大部分的原始信号，高频部分包含了大部分的噪声信号，文章中选择局部特性较好的DB3小波作为小波基函数，分解过程如下：
$$\begin{gathered} z^{\text{wdet }}={\downarrow }_2(z⊛\psi ) \\ {z}^{\text{wapp }}={\downarrow }_2(z⊛\phi ) \end{gathered}$$
其中：
ψ = ψ 1 ⊛ ψ 1 T φ = φ 1 ⊛ φ 1 T \psi=\psi_{1} \circledast\psi_{1}^{T}\\ \varphi=\varphi_{1} \circledast \varphi_{1}^{T} ψ=ψ1​⊛ψ1T​φ=φ1​⊛φ