運動模糊恢復
維納濾波—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 濾波後的均值，或者 拉普拉斯運算元 濾波的方差，或者取傅立葉變換後的高頻資訊均值。