運動模糊圖像復原研究的整體思路主要是用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像素時誤差就很大了，不能用作檢測
。。。。。。。。。。以上