运动模糊图像复原研究的整体思路主要是用matlab中的 imfilter（）函数对图像进行线性空间滤波，产生运动模糊图像，建立退化模型 → 通过radon变换来获取模糊参数，即点扩散函数PSF → 最后由估计得出的PSF再用维纳滤波对图像进行复原。由仿真实验得知，在已知PSF的情况下使用自相关函数的维纳滤波法对图像进行复原可以获得较好的复原效果，因此难点在于如何精确地估计运动模糊参数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 求出运动模糊方向。
实现代码：
1 close all; 2 clear all; 3 4 %% 读入并显示图像 5 filename = 'ex.jpg'; 6 I = imread(filename); 7 8 figure 9 imshow(uint8(I)); 10 title('原图'); 11 12 %% 生成运动模糊图像 13 PSF = fspecial('motion',50, 120); 14 g = imfilter(I, PSF, 'circular'); 15 figure 16 imshow(uint8(g)); 17 title('运动模糊图'); 18 19 %% 对运动模糊图像进行灰度化，并进行二维快速傅里叶变换，生成其频谱图 20 gb = rgb2gray(g); 21 figure 22 imshow(uint8(gb)); 23 PQ = paddedsize(size(gb)); 24 F = fft2(gb, PQ(1), PQ(2)); 25 figure 26 imshow(uint8(F)); 27 28 29 %% 将频谱压缩，居中 30 H = log(1+abs(F)); 31 Hc = fftshift(H); 32 figure 33 imshow(uint8(Hc)); 34 35 %% 用canny算子将压缩居中后的频谱图进行边缘检测，二值化 36 T = graythresh(Hc); 37 bw=edge(Hc, 'canny', T); 38 figure 39 imshow(bw); 40 41 %% 对二值化后的频谱图进行radon变换 42 theta = 1:180; 43 R = radon(bw, theta); 44 figure 45 imshow(R); 46 47 %% 计算出通过radon变换求出的模糊角度 48 MAX = max(max(R)); 49 [m, n] = find(R == MAX); 50 [M,N] = size(Hc); 51 beita = atan(tan(n*pi/180)*M/N)*180/pi;
（2）倒谱法：
1) 对模糊图像进行灰度化，并计算其二维傅里叶变换；
2) 对傅里叶变换取对数，平方后再进行一次反傅里叶变换得到其倒频谱；
3) 对倒频谱动态范围进行压缩；
4) 对压缩后的结果进行循环移位，使其低频成分居中；
5) 用canny算子对压缩居中后的频谱图像进行边缘检测使其二值化；
6) 将二值化后的频谱图做从1°~180°的radon变换；
7) 找出radon变换后的矩阵中的最大值，其对应的列数±90°即为所求的模糊方向。
实现代码：
1 close all; 2 clear all; 3 4 %% 读入并显示图像 5 filename = 'ex.jpg'; 6 I = imread(filename); 7 8 figure 9 imshow(uint8(I)); 10 title('原图'); 11 12 %% 生成运动模糊图像 13 PSF = fspecial('motion',80, 150); 14 g = imfilter(I, PSF, 'circular'); 15 figure 16 imshow(uint8(g)); 17 title('运动模糊图'); 18 19 %% 对运动模糊图像进行灰度化，并进行二维快速傅里叶变换，生成其频谱图 20 gb = rgb2gray(g); 21 figure 22 imshow(uint8(gb)); 23 PQ = paddedsize(size(gb)); 24 F = fft2(gb, PQ(1), PQ(2)); 25 figure 26 imshow(uint8(F)); 27 28 %% 作出倒频谱 29 F1 = log(1+abs(F)); 30 F2 = abs(F1).^2; 31 F3 = real(ifft2(F2)); 32 figure 33 imshow(uint8(F3)); 34 35 36 %% 将倒频谱压缩，居中 37 H = log(1+abs(F3)); % 将倒频谱动态范围进行压缩 38 Hc = fftshift(H); % 将压缩结果进行循环移位，使低频成分居中 39 figure 40 imshow(uint8(Hc)); 41 42 %% 通过阈值处理，边缘检测“canny”算子二值化倒频谱 43 T = graythresh(Hc); 44 bw=edge(Hc, 'canny', T); 45 figure 46 imshow(bw); 47 48 %% 对倒频谱从1°到180°作radon变换，以求出模糊角度 49 theta = 1:180; 50 R = radon(bw, theta); 51 figure 52 imshow(R); 53 54 %% 计算出通过倒频谱radon变换估计出的模糊角度 55 MAX = max(max(R)); 56 [m, n] = find(R == MAX); 57 if 90 < n <= 180 58 beita = n - 90; 59 else if 0 < n < 90 60 beita = n + 90; 61 else if n == [90;90] | n == [180;180] 62 beita = n(1); 63 end;
7、仿真结果及算法评价：
经过反复、多次、令人崩溃的调试参数，得出结论，模糊角度估计的精确度主要与边缘检测的阈值 T，运动模糊尺度和运动模糊方向这三个变量有关，经调试，阈值选取由graythresh 函数算出的全局阈值效果较好，而模糊尺度越小精度越低。
频域法：当模糊尺度大于20像素时，模糊方向处于0°~90°时，检测效果很好，误差低于0.5°，但当模糊尺度低于20像素时，由于中心化处理致使频谱图中间出现的十字亮线使模糊方向常为90°或180°；而当模糊尺度大于90°时，估计到的模糊方向像脱缰的野马。。完全不知道怎么才能估计出如此恶心的数据的。。
倒谱法：当模糊尺度大于40像素时，可以检测的模糊方向较大，基本从0°到180°都可以检测，不过误差在5°左右，模糊尺度小于40像素时误差就很大了，不能用作检测
。。。。。。。。。。以上