clear;clc
% 参考资料：http://blog.csdn.net/bluecol/article/details/47359421
% 参考数字图像处理(第三版) 冈萨雷斯 英文版公式(5.9-4)
% 约束最小二乘方滤波约束最小二乘方复原方法主要是针对有约束退化模型而言的，
% 这种图像复原方法只需有关噪声均值和方差的知识就能够对每一幅给定的图像进行
% 复原就以得到最优的效果。但有时该方程的解振荡的非常厉害。
I = im2double(imread('C:\Users\Administrator\Desktop\imread.png'));
tic
[m,n,~] = size(I);
subplot(1,3,1),imshow(I);
title('原图');
% Simulate a motion blur and additive noise.
LEN = 20;
THETA = 30;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');% 生成运动模糊图像
% psf2otf是把一个空间点扩散函数转换为频谱面的光学传递函数，其转换过程也是
% 一个二维的傅里叶变换。不过与fft2有区别
% psf2otf和fft2的区别：http://blog.csdn.net/bluecol/article/details/48288739
H = psf2otf(PSF,[m,n]);% 退化函数H(u,v)
noise_mean = 0;% 噪声均值
noise_var = 0.00001;% 噪声方差
blurred_noisy = imnoise(blurred,'gaussian',...
noise_mean, noise_var);% 加噪的运动模糊图像
subplot(1,3,2), imshow(blurred_noisy)
title('模糊退化且添加高斯噪声的图像')
% restoration using Home Made Constrained Least Squares Filtering.
p = [0 -1 0;-1 4 -1;0 -1 0];% 拉普拉斯算子
P = psf2otf(p,[m,n]);% 对拉普拉斯算子做傅立叶变换
gama = 0.001;% 公式里左边分式分母的一个参数
G = fft2(blurred_noisy);% 对模糊图像做傅立叶变换
numerator = conj(H);%公式里左边分式的分子，也就是H(u,v)的共轭
denominator = H.^2 + gama*(P.^2);%公式里左边分式的分母
deblurred2 = ifft2(numerator.*G./denominator);% 再从频域变换回时域,注意此处是点除
subplot(1,3,3), imshow(deblurred2)
title('约束最小二乘复原图');
toc
