保持结构不变的图像降噪

假定图像的区域是同构或者异构的。我们可以逐个处理每个像素，通过检测它的邻域结构类型（异构/同构）来估计像素的真实密度，从而减少噪声。

设 $X(p)$ 为真实图像, $Y(p)$ 为被高斯噪声污染的后的观测图像, $n(p)\in N(0,\sigma^2)$ 为均值为 0 ,方差为 $\sigma^2$ 的高斯噪声分布

X (p) = Y (p) + n (p)

$X(p)=Y(p)+n(p)$

一个同构的邻域 $N_{n\times n}(p)$ 假设是一个样本大小为 $N=n\times n$ ，均值为 $\mu$ ，方差为 $\sigma^2$ 的高斯随机变量。

一个异构的邻域则被假设为样本大小为N，满足两个高斯混合分布的变量，分布的先验为 $P_i$ ，均值为 $\mu_i$ ，以及相同的方差 $\sigma_i^2=\sigma^2$ , $i=1,2$ 。混合分布的概率密度函数为：

P (x) = \frac{1}{σ \sqrt{2 π}} \sum_{i = 0}^{1} P_{i} e x p {- \frac{(x - μ_{i})^{2}}{2 σ^{2}}}

$P(x)=\frac{1}{\sigma\sqrt{2\pi}}\sum_{i=0}^1{P_iexp{\{-\frac{(x-\mu_i)^2}{2\sigma^2}}}\}$

最大似然率（maximum likelihood radio test）为

如果 $\hat{S^2}\leq \sigma^2(1+C)$
邻域 $N_{n\times n}(p)$ 是同构的，其中 $\hat{S^2}$ 是邻域样本 $N_{n\times n}(p)$ 的方差：
$\hat{S^{2}} = \frac{1}{N} \sum_{q \in N_{n \times n} (p)} (Y (q) - \frac{1}{N} \sum_{r \in N_{n \times n (p)}} Y (r))^{2}$ $\hat{S^2}=\frac{1}{N}\sum_{q\in N_{n\times n}(p)}(Y(q)-\frac{1}{N}\sum_{r\in N_{n\times n(p)}}Y(r))^2$
参量C 由测试的显著性（sigificant level）决定 (例如，接受错误的同构的概率)。下表给出了C的值和对应的显著性。

如果邻域 $N_{n\times n}(p)$ 是同构的，则像素p的真实值的估计值为邻域所有像素的均值，此时的估计为高斯噪声下的最小无偏估计（最优估计）。即
$\hat{X} (p) = \frac{1}{N} \sum_{q \in N_{n \times n} (p)} Y (q)$ $\hat{X}(p)=\frac{1}{N}\sum_{q\in N_{n\times n}(p)}Y(q)$
如果邻域异构，则混合分布的先验和均值由前3个采样时刻估计：

\hat{μ_{i}} = \frac{1}{2} [β - (- 1)^{i} \sqrt{β^{2} - 4 λ}], \hat{P_{i}} = (- 1)^{i} \frac{\hat{u_{\bar{i}}} - c_{1}}{\hat{u_{\bar{i}}} - u_{i}}, i = 0, 1

$\hat{\mu_i}=\frac{1}{2}[\beta-(-1)^i\sqrt{\beta^2-4\lambda}], \hat{P_i}=(-1)^i\frac{\hat{u_{\bar{i}}}-c_1}{\hat{u_{\bar{i}}}-u_i},i=0,1$

\bar{i} = m o d (i + 1, 2) β = \frac{c_{3} - c_{1} c_{2}}{c_{2} - c_{1}^{2}}, λ = \frac{c_{1} c_{3} - c_{2}^{2}}{c_{2} - c_{1}^{2}}

$\bar{i}=mod(i+1,2) \beta=\frac{c_3-c_1c_2}{c_2-c_1^2},\lambda=\frac{c_1c_3-c_2^2}{c_2-c_1^2}$

c_{1} = m_{1}, c_{2} = m_{2} - σ^{2}, c_{3} = m_{3} - 3 m_{1} σ^{2}

$c_1=m_1,c_2=m_2-\sigma^2,c_3=m_3-3m_1\sigma^2$

$m_j,j=1,2,3$ 为 $N_{n\times n}(p)$ 的第j oder sample moments。

m_{j} = \frac{1}{N} \sum_{q \in N_{n \times n} (p)} Y^{j} (q)

$m_j=\frac{1}{N}\sum_{q\in N_{n\times n}(p)}Y^j(q)$
有了模型的估计参量，由贝叶斯方法可以得到两个混合分布的阈值

T_{c}

$T_c$

T_{c} = \frac{{\hat{μ}}_{0} + {\hat{μ}}_{1}}{2} + \frac{σ^{2}}{{\hat{μ}}_{1} - {\hat{μ}}_{0}} l n \frac{{\hat{P}}_{0}}{{\hat{P}}_{1}}

$T_c=\frac{\hat{\mu}_0+\hat{\mu}_1}{2}+\frac{\sigma^2}{\hat{\mu}_1-\hat{\mu}_0}ln\frac{\hat{P}_0}{\hat{P}_1}$
此时像素p的估计值为

{\hat{X}}_{p} = {\begin{cases} {\hat{μ}}_{1}, i f Y (p) > T_{c} \\ {\hat{μ}}_{0}, o t h e r w i s e \end{cases}

$\hat{X}_p=\begin{cases}\hat{\mu}_1, if Y(p)>T_c \cr\hat{\mu}_0, otherwise \end{cases}$
当信噪比

S N R = \frac{| μ_{0} - μ_{1} |}{σ}

$SNR=\frac{|\mu_0-\mu_1|}{\sigma}$ 大于2时，此算法能够有效的降噪同保持图像的结构。降噪的结果依赖于观测图像噪声方差

σ^{2}

$\sigma^2$ 的估计，以及 N, N可以定义或者利用噪声方差来估计。N越大，细的特征如线和角特征容易丢失。N越小，平滑效果越不明显。

matlab代码实例如下，完整代码(github)

RGB=imread('saturn.png');
I=rgb2gray(RGB);
noise=sqrt(0.025)*randn(size(I)) + 0;
NI=im2double(I)+noise;
figure,imshow(NI);title('stained by gaussian noise Image');
[SI,SNR]=ImageReduNoise(NI,0.025,'NeigborHoodSize',[7,7],'SigLevel','SigLevel3');
title('smooth Image')
K=wiener2(NI,[7,7]);
figure,imshow(K);title('wiener filter')

原图以及对原图加上方差为0.025的高斯噪声后的图像如下图所示:
original image
gaussian noise
当N取7，显著性为0.001时，降噪的结果如下图所示：
reduce noise
维纳自适应滤波器降噪结果如下图所示,可以看到效果要比维纳自适应滤波更好。
wiener filter

高斯噪声方差估计

此处讨论的是满足高斯分布的加性噪声
1. 利用均匀区域估计高斯噪声参数

选取图像中比较“平坦”的子带区域，此处被认为是均匀的，没有边缘存在，计算其方差和均值。将其直方图和估计出来的高斯分布形状比较。如果接近的，那么计算出的均值和方差就是高斯噪声的估计参数。

matlab代码实例如下，

I=imread('cameraman.tif');
% 加入高斯噪声方差为0.01，均值为0
bin=256;
noise=sqrt(0.01)*randn(size(I)) + 0;
% 图像转为double型和噪声叠加
NI=im2double(I)+noise;
figure,imshow(NI);
% 人工选择“平坦”区域
NIC=imcrop(NI);
% stastic mean and variance
[p,xvalue]=hist(double(NIC(:)),bin);
p=p./length(NIC(:));
u=sum(xvalue.*p);
var=sum(bsxfun(@minus,xvalue,u).^2.*p);
% 比较直方图和高斯分布
pdf=gaussmf(xvalue,[sqrt(var),u]).*((2*pi*var)^(-1/2));
figure,bar(pdf);
title('高斯分布');
figure,bar(p);title('直方图分布');

形状如图所示，两者很接近，此时高斯分布的方差估计值是0.0099和噪声实际方差0.01的误差为0.0001。均值估计值是0.0515，误差为0.0515。

基于滤波

我们知道边缘和噪声的二阶微分值都较大，但是不同的是噪声对laplacian算子比较敏感。因而可以用两个不同laplacian滤波器，对噪声图像滤波之后做差。由于边缘部分的滤波结果差不多，抵消了，剩下的就是噪声部分。两个Laplacian滤波器滤波后做差效果如下图所示：

自定义滤波器N滤波结果

具体算法：噪声图像imageI大小为[W,H]。两个Laplacian滤波器,L1和L2，估计噪声方差为 $\sigma_n$
L1L2

σ_{n} = \sqrt{\frac{π}{2}} \frac{1}{6 (W - 2) (H - 1)} \sum_{i m a g e I} | I * N |

$\sigma_n=\sqrt{\frac{\pi}{2}}\frac{1}{6(W-2)(H-1)} \sum_{imageI}|I*N|$
注意此方法只能估计均值为0的高斯分布噪声的方差。
matlab代码实例如下：

I=imread('cameraman.tif');
% 加入高斯噪声方差为0.01，均值为0
noise=sqrt(0.01)*randn(size(I)) + 0;
% 图像转为double型和噪声叠加
NI=im2double(I)+noise;
figure,subplot(121),imshow(I);title('原图');
subplot(122),imshow(NI);title('加噪声')
h1=fspecial('laplacian',0);%滤波器L1
h2=fspecial('laplacian',1);%滤波器L2
N=(h2-h1)*2;
imI=imfilter(NI,N);
%估计噪声方差
varn=sum(abs(imI(:)))./(size(imI,1)-2)./(size(imI,2)-2)./6.*sqrt(pi/2);
varn=varn.^2;

最后估计的varn值为0.0108，误差为0.0008

参考文献

[1]. Haris K, Efstratiadis S N, Maglaveras N, et al. Hybrid image segmentation using watersheds and fast region merging.[J]. IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society, 1998, 7(12):1684-1699.

[2]. Immerkaer J. Fast noise variance estimation[J]. Computer vision and image understanding, 1996, 64(2): 300-302.

[3]冈萨雷斯. 数字图像处理 (第二版)[J]. 电子工业出版社, 2007.

[4]Haris K, Tziritas G, Orphanoudakis S. Smoothing 2-D or 3-D images using local classification[C]//Proc. EUSIPCO. 1994