運動模糊圖像處理（一）----- 模糊角度估計的算法研究及matlab實現

時間 2019-12-13

標籤運動模糊圖像處理角度估計算法研究 matlab 實現欄目 MATLAB 简体版

原文原文鏈接

運動模糊圖像復原研究的總體思路主要是用matlab中的 imfilter（）函數對圖像進行線性空間濾波，產生運動模糊圖像，創建退化模型 → 經過radon變換來獲取模糊參數，即點擴散函數PSF → 最後由估計得出的PSF再用維納濾波對圖像進行復原。由仿真實驗得知，在已知PSF的狀況下使用自相關函數的維納濾波法對圖像進行復原能夠得到較好的復原效果，所以難點在於如何精確地估計運動模糊參數PSF。算法

一、基本原理：函數

點擴散函數PSF主要有兩個重要參數：（1）模糊方向；（2）模糊尺度。本次主要是針對第一個參數----模糊方向的估計進行了研究。運動模糊方向是指運動方向與水平方向的夾角，由文獻得知運動模糊主要是下降了運動方向的高頻成分，而對其餘方向的高頻成分影響較小。常見的辨識方法有頻域法和倒譜法，wym 兩種方法都試過，仿真實驗結果表兩種方法各有好處。ui

頻域法的原理是將退化圖像進行二維傅里葉變換，獲得具備相互平行的規則明暗條紋的頻譜。設暗紋與 x 軸正向夾角爲 φ ，運動模糊方向與 x 軸夾角爲 θ ，圖像尺寸爲 M × N,根據傅里葉變換的時頻特性能夠知道，可經過公式 tan(θ) = tan(φ − 90°) × M/N 獲得模糊角度 θ ,所以只要經過 Radon 變換檢測出頻譜暗條紋與水平方向的夾角便可到運動模糊方向。spa

倒譜法的主要原理是先將退化圖像進行二維傅里葉變換，而後取對數，再進行反傅里葉變換獲得退化圖像的倒頻譜,分離出退化圖像的模糊信息，進而經過 Radon 變換獲得運動模糊方向。調試

Radon 變換是對頻譜圖上某一指定角度進行線積分，經過計算1°~180°的Radon變換獲得180列的矩陣 R，每一列向量是圖像在一個角度上沿一族直線的積分投影，由於積分直線束與頻譜中的亮暗條紋平行，因此所得的投影向量中應有一個最大值，在頻域法中最大值所對應的列數就等於模糊方向與x軸正方向水平夾角；在倒譜法中，最大值對應的列數 ±90°即爲所求的模糊角度。code

具體理論和公式推導就不列出來了。。有興趣的同窗請 STFW。。（什麼？不知道什麼是 STFW? 請自行 STFW。。）blog

二、算法實現：ci

（1）頻域法it

1) 對模糊圖像進行灰度化，並計算其二維傅里葉變換；io

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;