维纳滤波—deconvwnr函数利用维纳滤波器来对图像模糊修复
function image_restoration_deconvwnr()
%Read image
I = im2double(imread('lena.tif'));
I=rgb2gray(I);
figure,subplot(2,3,1),imshow(I);
title('Original Image');
%Simulate a motion blur
LEN = 21;
THETA = 11;
PSF = fspecial('motion', LEN, THETA); %获得滤波算子 PSF
blurred = imfilter(I, PSF, 'conv', 'circular');
subplot(2,3,2),imshow(blurred);
title('Blurred Image');
%Restore the blurred image
wnr1 = deconvwnr(blurred, PSF, 0);
subplot(2,3,3),imshow(wnr1);
title('Restored Image');
%Simulate blur and noise
noise_mean = 0;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, 'gaussian', ...
noise_mean, noise_var);
subplot(2,3,4),imshow(blurred_noisy)
title('Simulate Blur and Noise')
%Restore the blurred and noisy image:First attempt
wnr2 = deconvwnr(blurred_noisy, PSF, 0);
subplot(2,3,5);imshow(wnr2);title('Restoration of Blurred, Noisy Image Using NSR = 0')
%Restore the Blurred and Noisy Image: Second Attempt
signal_var = var(I(:));
wnr3 = deconvwnr(blurred_noisy, PSF, noise_var / signal_var);
subplot(2,3,6),imshow(wnr3)
title('Restoration of Blurred, Noisy Image Using Estimated NSR');
end维纳滤波需要估计图像的信噪比(SNR)或者噪信比(NSR)，信号的功率谱使用图像的方差近似估计，噪声分布是已知的。从第一排中可以看出，若无噪声，此时维纳滤波相当于逆滤波，恢复运动模糊效果是极好的。从第二排可以看出噪信比估计的准确性对图像影响比较大的，二排中间效果几乎不可用。
逆滤波
function image_restoration_matlab()
% Display the original image.
I = im2double(imread('lena.tif'));
I=rgb2gray(I);
[hei,wid,~] = size(I);
subplot(2,3,1),imshow(I);
title('Original Image (courtesy of MIT)');
% Simulate a motion blur.
LEN = 21;
THETA = 11;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');
subplot(2,3,2), imshow(blurred); title('Blurred Image');
% Inverse filter
If = fft2(blurred);
Pf = fft2(PSF,hei,wid);
deblurred = ifft2(If./Pf);
subplot(2,3,3), imshow(deblurred); title('Restore Image')
% Simulate additive noise.
noise_mean = 0;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, 'gaussian', ...
noise_mean, noise_var);
subplot(2,3,4), imshow(blurred_noisy)
title('Simulate Blur and Noise')
% Try restoration assuming no noise.
If = fft2(blurred_noisy);
deblurred2 = ifft2(If./Pf);
subplot(2,3,5), imshow(deblurred2)
title('Restoration of Blurred Assuming No Noise');
% Try restoration with noise is known.
noisy = blurred_noisy - blurred;
Nf = fft2(noisy);
deblurred2 = ifft2(If./Pf - Nf./Pf);
subplot(2,3,6), imshow(deblurred2)
title('Restoration of Blurred with Noise Is Known')
end逆滤波对噪声非常敏感，除非我们知道噪声的分布情况（事实上，这也很难知道），逆滤波几乎不可用，可以从二排中间看出，恢复图像效果极差。但若知道噪声分布，也是可以完全复原信息的。可以从二排最后一张图可以看出。
注：图像恢复时需要用到 PSF 卷积核,该卷积核为生成运动模糊得卷积核，在实际应用场景是不知道的。
如何求解： PSF
PSF 是由 LEN 和 Theta 决定的，如果能求出 运动模糊的方向，和对应的偏移量就能求出 PSF
1、基本原理：
点扩散函数PSF主要有两个重要参数：（1）模糊方向；（2）模糊尺度。本次主要是针对第一个参数----模糊方向的估计进行了研究。运动模糊方向是指运动方向与水平方向的夹角，由文献得知运动模糊主要是降低了运动方向的高频成分，而对其他方向的高频成分影响较小。常见的辨识方法有频域法和倒谱法，wym 两种方法都试过，仿真实验结果表两种方法各有好处。
频域法的原理是将退化图像进行二维傅里叶变换，得到具有相互平行的规则明暗条纹的频谱。设暗纹与 x 轴正向夹角为 φ ，运动模糊方向与 x 轴夹角为 θ ，图像尺寸为 M × N,根据傅里叶变换的时频特性可以知道，可通过公式 tan(θ) = tan(φ − 90°) × M/N 得到模糊角度 θ ,因此只要通过 Radon 变换检测出频谱暗条纹与水平方向的夹角即可到运动模糊方向。
倒谱法的主要原理是先将退化图像进行二维傅里叶变换，然后取对数，再进行反傅里叶变换得到退化图像的倒频谱,分离出退化图像的模糊信息，进而通过 Radon 变换得到运动模糊方向。
Radon 变换是对频谱图上某一指定角度进行线积分，通过计算1°~180°的Radon变换得到180列的矩阵 R，每一列向量是图像在一个角度上沿一族直线的积分投影，因为积分直线束与频谱中的亮暗条纹平行，所以所得的投影向量中应有一个最大值，在频域法中最大值所对应的列数就等于模糊方向与x轴正方向水平夹角；在倒谱法中，最大值对应的列数 ±90°即为所求的模糊角度。
具体理论和公式推导就不列出来了。。有兴趣的同学请 STFW。。（什么？不知道什么是 STFW? 请自行 STFW。。）
2、算法实现：
（1）频域法
1) 对模糊图像进行灰度化，并计算其二维傅里叶变换；
2) 对傅里叶变换值的动态范围进行压缩；
3) 对压缩后的结果进行循环移位，使其低频成分居中；
4) 用canny算子对压缩居中后的频谱图像进行边缘检测使其二值化；
5) 将二值化后的频谱图做从1°~180°的radon变换；
6) 找出radon变换后的矩阵中的最大值，求出其对应的列数 n；
7) 通过公式 tan(θ) = tan(φ − 90°) × M/N = tan(n − 90°) × M/N 求出运动模糊方向。
实现代码：
close all;
clear all;
%% 读入并显示图像
filename = 'ex.jpg';
I = imread(filename);
figure
imshow(uint8(I));
title('原图');
%% 生成运动模糊图像
PSF = fspecial('motion',50, 120);
g = imfilter(I, PSF, 'circular');
figure
imshow(uint8(g));
title('运动模糊图');
%% 对运动模糊图像进行灰度化，并进行二维快速傅里叶变换，生成其频谱图
gb = rgb2gray(g);
figure
imshow(uint8(gb));
PQ = paddedsize(size(gb));
F = fft2(gb, PQ(1), PQ(2));
figure
imshow(uint8(F));
%% 将频谱压缩，居中
H = log(1+abs(F));
Hc = fftshift(H);
figure
imshow(uint8(Hc));
%% 用canny算子将压缩居中后的频谱图进行边缘检测，二值化
T = graythresh(Hc);
bw=edge(Hc, 'canny', T);
figure
imshow(bw);
%% 对二值化后的频谱图进行radon变换
theta = 1:180;
R = radon(bw, theta);
figure
imshow(R);
%% 计算出通过radon变换求出的模糊角度
MAX = max(max(R));
[m, n] = find(R == MAX);
[M,N] = size(Hc);
beita = atan(tan(n*pi/180)*M/N)*180/pi;（2）倒谱法：
1) 对模糊图像进行灰度化，并计算其二维傅里叶变换；
2) 对傅里叶变换取对数，平方后再进行一次反傅里叶变换得到其倒频谱；
3) 对倒频谱动态范围进行压缩；
4) 对压缩后的结果进行循环移位，使其低频成分居中；
5) 用canny算子对压缩居中后的频谱图像进行边缘检测使其二值化；
6) 将二值化后的频谱图做从1°~180°的radon变换；
7) 找出radon变换后的矩阵中的最大值，其对应的列数±90°即为所求的模糊方向。
实现代码：
close all;
clear all;
%% 读入并显示图像
filename = 'ex.jpg';
I = imread(filename);
figure
imshow(uint8(I));
title('原图');
%% 生成运动模糊图像
PSF = fspecial('motion',80, 150);
g = imfilter(I, PSF, 'circular');
figure
imshow(uint8(g));
title('运动模糊图');
%% 对运动模糊图像进行灰度化，并进行二维快速傅里叶变换，生成其频谱图
gb = rgb2gray(g);
figure
imshow(uint8(gb));
PQ = paddedsize(size(gb));
F = fft2(gb, PQ(1), PQ(2));
figure
imshow(uint8(F));
%% 作出倒频谱
F1 = log(1+abs(F));
F2 = abs(F1).^2;
F3 = real(ifft2(F2));
figure
imshow(uint8(F3));
%% 将倒频谱压缩，居中
H = log(1+abs(F3)); % 将倒频谱动态范围进行压缩
Hc = fftshift(H); % 将压缩结果进行循环移位，使低频成分居中
figure
imshow(uint8(Hc));
%% 通过阈值处理，边缘检测“canny”算子二值化倒频谱
T = graythresh(Hc);
bw=edge(Hc, 'canny', T);
figure
imshow(bw);
%% 对倒频谱从1°到180°作radon变换，以求出模糊角度
theta = 1:180;
R = radon(bw, theta);
figure
imshow(R);
%% 计算出通过倒频谱radon变换估计出的模糊角度
MAX = max(max(R));
[m, n] = find(R == MAX);
if 90 < n <= 180
beita = n - 90;
else if 0 < n < 90
beita = n + 90;
else if n == [90;90] | n == [180;180]
beita = n(1);
end;7、仿真结果及算法评价：
经过反复、多次、令人崩溃的调试参数，得出结论，模糊角度估计的精确度主要与边缘检测的阈值 T，运动模糊尺度和运动模糊方向这三个变量有关，经调试，阈值选取由graythresh 函数算出的全局阈值效果较好，而模糊尺度越小精度越低。
频域法：当模糊尺度大于20像素时，模糊方向处于0°~90°时，检测效果很好，误差低于0.5°，但当模糊尺度低于20像素时，由于中心化处理致使频谱图中间出现的十字亮线使模糊方向常为90°或180°；而当模糊尺度大于90°时，估计到的模糊方向像脱缰的野马。。完全不知道怎么才能估计出如此恶心的数据的。。
倒谱法：当模糊尺度大于40像素时，可以检测的模糊方向较大，基本从0°到180°都可以检测，不过误差在5°左右，模糊尺度小于40像素时误差就很大了，不能用作检测
关于偏移距离的计算 可以使用 折中法试探，每次试探一半的距离，然后根据恢复结果的图像质量进一步选取新的距离值。 关于图像质量评价可以使用纹理特征，例如 Sobel 滤波后的均值，或者 拉普拉斯算子 滤波的方差，或者取傅里叶变换后的高频信息均值。
