改進維納濾波的實作——光學稀疏孔徑成像系統圖像恢複算法研究 陳灏

        本文簡單實作了浙江大學陳灏碩士學位論文《光學稀疏孔徑成像系統圖像恢複算法研究》(2017)中4.2章節提及的改進維納濾波算法。

        前情提要：維納濾波函數中不同K(噪聲與原屬圖像功率譜之比)的取值對應複原效果不同，當K值比較小時(10^-4~10^-3)時，圖像邊緣存在着嚴重的振鈴效應，而随着K的慢慢增大，振鈴效應越來越弱直至消失，最後圖像變得越來越模糊。是以作者提出了一中基于PSNR的圖像品質評價體系的改進的維納濾波算法，目的是選取最優的K值，使得用該K做維納濾波得到的複原後的圖像與原始圖像最為接近（PSNR值最大）。

        原文在正文第36頁：

        但原文似乎沒有對此方法的詳細介紹，是以姑且用一個簡單的方法考慮：每次K按一定方式，從10^-5~到10^3量級增加，判斷以目前K值進行維納濾波得到的複原圖像與原始圖像的差異度（即PSNR）是否比曆史最優值更好，類似于周遊找最大值的政策，但量級跳度過大，如果按照 10^-5~到10^3量級的原則來的話，需要進行10^8次計算，這無疑增加了計算量。根據觀察，可以根據原文中K值與PSNR值曲線來進行K的跳度限制：最開始随着K值的增加，PSNR迅速增大到最大值後開始越來越平穩地減小，相對應可設定最開始階段K值的增加（步長）為10^-5級，之後根據K值來對步長進行擴增，如10倍等。當然精度可以根據需求或者個人喜好等随意修改，在代碼中，初始化K為0，步長step為1*10^-5（matlab中可以簡單表示為1e-5），當K值增加到是步長的100倍時，步長相應地增加到原來的10倍，直至K大于1000，停止循環。

      matlab中，可直接使用函數deconvwnr進行維納濾波操作。為友善起見，将同系函數的使用方法進行總結：

調用方式	方法名	參數解釋
J = deconvwnr(I,psf) J = deconvwnr(I,psf,nsr) J = deconvwnr(I,psf,ncorr,icorr)	維納濾波	psf:點擴充函數(卷積核) nsr:加性噪聲的噪信比(K) ncorr:噪聲自相關函數 icorr:I的自相關函數
J = deconvlucy(I,psf) J = deconvlucy(I,psf,iter) J = deconvlucy(I,psf,iter,dampar) J = deconvlucy(I,psf,iter,dampar,weight) J = deconvlucy(I,psf,iter,dampar,weight,readout) J = deconvlucy(I,psf,iter,dampar,weight,readout,subsample)	Lucy-Richardson濾波	iter:疊代次數 dampar:阻尼門檻值  weight:權重矩陣 readout:讀出相機噪聲 subsample:子采樣倍數
[J,psfr] = deconvblind(I,psfi) [J,psfr] = deconvblind(I,psfi,iter) [J,psfr] = deconvblind(I,psfi,iter,dampar) [J,psfr] = deconvblind(I,psfi,iter,dampar,weight) [J,psfr] = deconvblind(I,psfi,iter,dampar,weight,readout) [J,psfr] = deconvblind(___,fun)	盲去卷積	psfi：點擴充函數的初始估計 psfr:複原後的點擴充函數 fun：描述點擴充函數附加限制的函數句柄
J = deconvreg(I,psf) J = deconvreg(I,psf,np) J = deconvreg(I,psf,np,lrange) J = deconvreg(I,psf,np,lrange,regop) [J,lagra] = deconvreg(___)	最小二乘濾波	np：加性噪聲強度 lrange：搜尋最優解範圍 regop：正則算子 lagra：拉格朗日乘子

        psf相當于卷積核，由檔案讀取，也可以通過其他方式設定。psnr可以直接調用matlab函數，亦可以手動實作。data為n*4維矩陣，每一行存儲了每次計算對應的次數、K值、步長和PSNR值。

        代碼和注釋如下，相應條件可自行調整：

clc
clear
close all
f=rgb2gray(imread('timg.jpg'));
[m,n]=size(f);
Psf=load('psfdata.txt');
%   C=uint8(conv2(f,Psf,'same')); %另一種卷積方式
C = imfilter(f,Psf,'conv','circular'); %卷積得到模糊圖像
K=0;
step=1e-5; %步長 
Img=[];%最大PSNR對應的圖像
best_PNSR=-Inf;%最大值
best_K=-Inf;%最大值對應的K
t=1;%計算次數
while K<=1000
    img=deconvwnr(C,Psf,K);%維納濾波
    new_PSNR=psnr(img,f);%計算濾波複原後的圖像與原圖的PSNR值
    data(t,1:4)=[t,K,step,new_PSNR];%儲存次數、K值、步長、PSNR值
    if new_PSNR>best_PNSR%最大PSNR對應的更新
        best_PNSR=new_PSNR;
        best_K=K;
        Img=img;
    end
    if K/step>100%步長的更新
        step=step*10;
    end
    K=K+step;%K值的更新
    t=t+1;
end
%對K值和PSNR值進行繪制，由于圖像前後差距過大，隻取前120行
x=data(1:120,2);
y=data(1:120,4);
plot(x,y)
xlabel('K值')
ylabel('PSNR')
figure,
subplot(131),imshow(f),title('原始圖像')
subplot(132),imshow(C),title('模糊圖像')
subplot(133),imshow(Img),title('複原圖像')
           

執行結果如下： 

模糊圖像PSNR:  19.9984

複原圖像PSNR:  27.3320