保持结构不变的图像降噪及高斯噪声估计

github地址

保持结构不变的图像降噪

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

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

$$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}{\sigma\sqrt{2\pi}}\sum_{i=0}^1{P_iexp{\{-\frac{(x-\mu_i)^2}{2\sigma^2}}}\}$$

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

1. 如果 $\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$$
   参量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)$$
2. 如果邻域异构，则混合分布的先验和均值由前3个采样时刻估计：

$$\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}=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-\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)$$
有了模型的估计参量，由贝叶斯方法可以得到两个混合分布的阈值 $T_c$
$$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{\mu}_1, if Y(p)>T_c \cr\hat{\mu}_0, otherwise \end{cases}$$
当信噪比 $SNR=\frac{|\mu_0-\mu_1|}{\sigma}$大于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的高斯噪声后的图像如下图所示:


当N取7，显著性为0.001时，降噪的结果如下图所示：

维纳自适应滤波器降噪结果如下图所示,可以看到效果要比维纳自适应滤波更好。

高斯噪声方差估计

此处讨论的是满足高斯分布的加性噪声
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。

1. 基于滤波

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

具体算法：噪声图像imageI大小为[W,H]。两个Laplacian滤波器,L1和L2，估计噪声方差为$\sigma_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