图像复原之约束最小二乘方滤波

图像复原，简单讲，就是恢复图像本来的面貌，但由于各种原因如图像采集过程中出现的误差导致得到的数字图像不清晰，不是我们人眼看到的实物场景那样，因此需要采取技术手段去除图像的不清晰。约束最小二乘方滤波就是其中一种较好的方法。在维纳滤波那一篇讲过，维纳滤波要求未退化图像和噪声的功率谱必须是已知的，通常这两个功率谱很难估计，尽管用一个常数去估计功率谱比，然而并不总是一个合适的解。约束最小二乘方滤波要求噪声的方差和均值，这些参数可通过给定的退化图像计算出来，这是约束最小二乘方滤波的一个重要优点。
 一个图像采集系统输入输出的关系可以表示为g(x,y)=H[f(x,y)]+η(x,y)，表达为向量-矩阵形式为g=Hf+η，明确地以矩阵形式来表达问题可以简化复原技术的推导。约束最小二乘方滤波的核心是解决退化函数H对噪声的敏感性问题，而减少噪声敏感性问题的一种方法是以平滑度量的最佳复原为基础的，如图像的二阶导数即拉普拉斯变换。于是，我们找到一个带约束条件的最小准则函数C，定义如下：

其约束条件为

其中是欧几里德向量范数，是未退化图像的估计。是拉普拉斯算子。求函数C的最小值，便得到最好的平滑效果即最佳复原。怎样找到最小值呢，这里我们采用拉格朗日乘数法，在频率域中，函数C可表示为C=，其中P是拉普拉斯算子的傅里叶变换。则频率域中，拉格朗日函数为

其中N是加性噪声η的傅里叶变换，对求导，得到的最小值表达式如下：

关于gamma值的确定：

令r = g-Hf ，(r)=

po出结论，(r)是gamma的单调递增函数，于是乎，在我们有了退化函数，退化后的图像，我们可以令

我们通过调整gamma的值令上式近乎相等，那么得到的gamma就是基本模拟了噪声的真实分布，对噪声的估计，我们通过噪声的均值和方差就可以得到了：

为什么可以这么估计，一个式子就可以解释了：

于是乎，我们只通过已知噪声的均值与方差，我们就可以实现对原始图像的恢复，而不需要像维纳函数那样还必须知道噪声的功率谱，得到噪声的真实分布在实际中并不是很现实。

matlab代码：

f = imread('test.tif');
[h,w]=size(f);
subplot(2,3,1),imshow(f);
title('原始图像')
PSF = fspecial('motion',21,11);%点扩散函数，运动模糊
blurred = imfilter(f,PSF,'conv','circular');
subplot(2,3,2),imshow(blur_f);
title('运动模糊的图像');
noise_mean=0;
noise_var = 0.0001;
g = imnoise(blurred, 'gaussian',noise_mean, noise_var);
subplot(2,3,3),imshow(g);
title('加了噪声的模糊图像');
H = psf2otf(PSF,[h,w]);
p = [0 -1 0;-1 4 -1;0 -1 0];%拉普拉斯模板
P = psf2otf(p,[h,w]);
G = fft2(g);
gamma=0.01;
numerator = conj(H)
denominator = H.^2+gamma*(P.^2);
deblur_f = ifft2(numerator.*G./denominator);
subplot(2,3,4);
imshow(deblur_f,[]);
title('约束最小二乘复原的图像')
%迭代计算最优的gamma
while ((g-blurred).^2-h*w*(noise_mean.^2+noise_mean.^2))>0.00001
    gamma=gamma-0.001
end
numerator = conj(H)
denominator = H.^2+gamma*(P.^2);
deblur_f = ifft2(numerator.*G./denominator);
subplot(2,3,5);
imshow(deblur_f,[]);
title('约束最小二乘复原最佳gamma的图像')