ch8-圖象復(fù)原2

1、 逆濾波比較簡(jiǎn)單，但沒(méi)有清楚地說(shuō)明如何處逆濾波比較簡(jiǎn)單，但沒(méi)有清楚地說(shuō)明如何處理噪聲，而維納濾波綜合了退化函數(shù)和噪聲理噪聲，而維納濾波綜合了退化函數(shù)和噪聲統(tǒng)計(jì)特性兩個(gè)方面進(jìn)行復(fù)原處理統(tǒng)計(jì)特性兩個(gè)方面進(jìn)行復(fù)原處理.5 維納濾波器維納濾波器逆濾波方法不能完全恢復(fù)原始信號(hào)逆濾波方法不能完全恢復(fù)原始信號(hào)f(x,y)，而只能，而只能求出求出f(x,y)的一個(gè)估計(jì)值的一個(gè)估計(jì)值f x,y希望找到一種方法，在有噪聲條件下，從退化圖像希望找到一種方法，在有噪聲條件下，從退化圖像g(x,y)復(fù)原出復(fù)原出f(x,y)的估計(jì)值，該估計(jì)值符合一定的準(zhǔn)則。的估計(jì)值，該估計(jì)值符合一定的準(zhǔn)則。維納濾波維納濾波 (Wiene

2、r filtering)=最小均方差濾波最小均方差濾波l維納濾波是最常用的圖像恢復(fù)方法維納濾波是最常用的圖像恢復(fù)方法l基于維納濾波的圖像恢復(fù)方法是基于維納濾波的圖像恢復(fù)方法是19671967年提出的年提出的C.W. Helstrom, “Image restoration by the method of lest sqaures,” Journal of the Optical Scoiety of America, vol.57, no.3, pp.297-303, 1967.C.W.Helstrom, This weeks citation classic, 19821967-1982年

3、年SCI引用超過(guò)引用超過(guò)125次次.N.Wiener, “The extrapolation, interpolation and smoothing of stationary time series”, New York: Wiely, 1949.22min: eEf x,yf x,y在均方誤差值最小的準(zhǔn)則下得到的在均方誤差值最小的準(zhǔn)則下得到的 稱為對(duì)稱為對(duì)f(x,y)的最小二乘方估計(jì)。的最小二乘方估計(jì)。f x,y按照該準(zhǔn)則得到的濾波器叫維納濾波器。按照該準(zhǔn)則得到的濾波器叫維納濾波器。目標(biāo)：目標(biāo)：使得復(fù)原后圖像使得復(fù)原后圖像 與原始圖像與原始圖像 的均方的均方誤差最?。赫`差最?。?,(yx

4、ff x,y因此維納濾波器又稱為最小均方差濾波器因此維納濾波器又稱為最小均方差濾波器線性濾波：尋找點(diǎn)擴(kuò)散函數(shù)hw(x,y)，使得),(*),(),(yxgyxhyxfw),(),(),(vuGvuHvuFW則有則有由由Andrews和和Hunt推導(dǎo)滿足這一要求的傳遞函數(shù)為：推導(dǎo)滿足這一要求的傳遞函數(shù)為： 2*( , )( , )( , )( , )( , )wnfHu vHu vS u vH u vSu v2*( , )( , )( , )( , )( , )/( , )nfHu vF u vG u vH u vS u vSu v這里，這里， 是成像系統(tǒng)傳遞函數(shù)是成像系統(tǒng)傳遞函數(shù)H(u,vH(

5、u,v) )的復(fù)共軛；的復(fù)共軛；S Sn( (u,v) 是噪聲功率譜；是噪聲功率譜；S Sf (u,v)是輸入圖像的功率譜。是輸入圖像的功率譜。 2*( , )( , )( , )( , )( , )/( , )nfHu vF u vG u vH u vS u vSu v221|( , )|( , )( , )( , ) |( , )|( , )/( , )nfH u vF u vG u vH u vH u vS u vSu v或者：*( , )Hu v2( , )( , )nS u vN u v2( , )( , )fSu vF u v 計(jì)算退化圖像計(jì)算退化圖像g(x,y)的二維的二維Fou

6、rier變換變換G(u,v)Wiener濾波的過(guò)程濾波的過(guò)程: 計(jì)算理想圖像的頻譜估計(jì)計(jì)算理想圖像的頻譜估計(jì) 計(jì)算點(diǎn)擴(kuò)展函數(shù)計(jì)算點(diǎn)擴(kuò)展函數(shù)h(x,y)的二維的二維Fourier變換變換H(u,v) 計(jì)算退化計(jì)算退化 圖像和噪聲的功率譜圖像和噪聲的功率譜Sf(u,v),Sn(u,v)( , )( , ) ( , )WF u vHu v G u v 計(jì)算濾波器計(jì)算濾波器HW(u,v) 求反求反Fourier變換變換2*( , )( , )( , )( , )( , )wnfHu vHu vS u vH u vSu v這一方法有如下特點(diǎn)：這一方法有如下特點(diǎn)：（1）當(dāng)）當(dāng)H（u,v）0或幅值很小時(shí)，分

7、母不為零，不會(huì)造成嚴(yán)重或幅值很小時(shí)，分母不為零，不會(huì)造成嚴(yán)重的運(yùn)算誤差。的運(yùn)算誤差。（2）在信噪比高的圖像中，即）在信噪比高的圖像中，即Sn(u,v)Sf(u,v)維納濾波復(fù)原法特點(diǎn)：1(,)(,)WHuvHuv2*( , )( , )( , )( , )( , )/( , )nfHu vF u vG u vH u vS u vSu v如果沒(méi)有噪聲，就成為逆濾波如果沒(méi)有噪聲，就成為逆濾波(3) 當(dāng)理想圖像功率譜當(dāng)理想圖像功率譜Sf (u,v)0)時(shí)時(shí) ,表明我們不可能表明我們不可能從全是噪聲的圖像中恢復(fù)出任何有意義的信號(hào)從全是噪聲的圖像中恢復(fù)出任何有意義的信號(hào) ( , )0F u v（4）往往

8、未退化圖像的功率譜）往往未退化圖像的功率譜Sf (u,v)難以知道，用下式近似表示：難以知道，用下式近似表示：),(),(),(),(1),(22vuGKvuHvuHvuHvuF維納濾波函數(shù)：測(cè)試維納濾波效果：逆濾波和維納濾波的比較v維納濾波的結(jié)果非常接近原始圖像,比逆濾波要好全逆濾波的結(jié)果半徑受限的逆濾波結(jié)果維納濾波的結(jié)果 (交互選擇K)逆濾波和維納濾波的比較v(a)運(yùn)動(dòng)模糊及均值為0方差為650的加性高斯噪聲污染的圖像 v(b) 逆濾波的結(jié)果v(c) 維納濾波的結(jié)果v(d)-(f) 噪聲幅度的方差比(a)小1個(gè)數(shù)量級(jí)v(g)-(i) 噪聲幅度的方差比(a)小5個(gè)數(shù)量級(jí)逆濾波逆濾波(inve

9、rse filtering)線性代數(shù)恢復(fù)線性代數(shù)恢復(fù)圖像恢復(fù)方法圖像恢復(fù)方法頻域圖像恢復(fù)頻域圖像恢復(fù)維納濾波維納濾波(Wiener filtering)代數(shù)逆濾波代數(shù)逆濾波 奇異值偽逆濾波。奇異值偽逆濾波。線性代數(shù)恢復(fù)在線性代數(shù)恢復(fù)在1970s1970s提出提出B.R.Hunt, “A matrix theory proof the discrete convolution theorem,” IEEE Trans. Audio Electroacoust , vol.19, no.4, pp.285-288, 1971.B.R.Hunt, “The application of constr

10、ained least squares estimation to image restoration by digital computer,” IEEE Trans. Computers, vol.22, no.9, pp.805-812, 1973. 約束最小二乘濾波約束最小二乘濾波gHfn 投影法。投影法。6 代數(shù)逆濾波代數(shù)逆濾波無(wú)約束代數(shù)逆濾波無(wú)約束代數(shù)逆濾波 (unconstrained restoration)已知退化圖像的向量形式已知退化圖像的向量形式g和退化矩陣和退化矩陣H, ,則無(wú)約則無(wú)約束逆濾波恢復(fù)的的圖像為束逆濾波恢復(fù)的的圖像為: :if is a square mat

11、rix11TTHfH HH gH g結(jié)果結(jié)果設(shè)恢復(fù)的圖像為設(shè)恢復(fù)的圖像為 ，如果不考慮噪聲，用，如果不考慮噪聲，用 表示恢表示恢復(fù)誤差。則代數(shù)逆濾波的目的是最小化目標(biāo)函數(shù)復(fù)誤差。則代數(shù)逆濾波的目的是最小化目標(biāo)函數(shù)f( )efgHf22( ) | ( )|() ()TWeffgHfgHfgHf2argmin |ffgHf即即推導(dǎo)推導(dǎo):( )2()0 fHgHffTW使導(dǎo)數(shù)為零使導(dǎo)數(shù)為零H gH HfTT0H gH HfTT1fH HH gTT無(wú)約束代數(shù)逆濾波無(wú)約束代數(shù)逆濾波 (unconstrained restoration)2argmin( , )( , )F u vF u vdudvFF注

12、意注意Wiener濾波的差別濾波的差別約束最小二乘濾波約束最小二乘濾波 (constrained least squares restoration)22|gHfn相對(duì)于無(wú)約束最小二乘濾波，約束最小二乘濾在最小化目標(biāo)函相對(duì)于無(wú)約束最小二乘濾波，約束最小二乘濾在最小化目標(biāo)函數(shù)時(shí)，利用了估計(jì)出的噪聲的信息作為約束條件。而我們已經(jīng)數(shù)時(shí)，利用了估計(jì)出的噪聲的信息作為約束條件。而我們已經(jīng)知道噪聲模型等信息可以較容易地從圖像中估計(jì)出來(lái)。通常把知道噪聲模型等信息可以較容易地從圖像中估計(jì)出來(lái)。通常把噪聲的范數(shù)作為約束，即噪聲的范數(shù)作為約束，即基本思想基本思想22| =| |gHfn令令Q為為f的線性矩陣算子，

13、約束最小二乘恢復(fù)問(wèn)題就是在滿足約的線性矩陣算子，約束最小二乘恢復(fù)問(wèn)題就是在滿足約束條件束條件 下，使下，使 最小化的問(wèn)題。最小化的問(wèn)題。2|Qf22 |s.t| =|. |gHfn2 |min|Qf1fH HQ QH gTTT例如選擇例如選擇Q為有限差分矩陣，使得二階差分的能量最小為有限差分矩陣，使得二階差分的能量最小211200( , )MNxyPf x y約束最小二乘濾波約束最小二乘濾波 (constrained least squares restoration)f (x, y)在在(x, y)處的二階微分處的二階微分)1,() 1,(), 1(), 1(),(4),(22222yxfy

14、xfyxfyxfyxfyfxfyxf222fQgHfnminL( )f 12201TTTTL( f )Q QfH ( gHf )ff( H HQ Q ) H g 其其中中用拉格朗日法求用拉格朗日法求微分，微分，可以用來(lái)調(diào)節(jié)以滿足約束條件?？梢杂脕?lái)調(diào)節(jié)以滿足約束條件。1fH HQ QH gTTT約束最小二乘濾波約束最小二乘濾波 (constrained least squares restoration)22 |s.t| =|. |gHfn2 |min|Qf推導(dǎo)推導(dǎo)clearF=checkerboard(8);subplot(2,2,1),imshow(F,);title(a) 原始圖像,Fon

15、tSize,12)PSF=fspecial(motion,7,45);MF=imfilter(F,PSF,circular);noise=imnoise(zeros(size(F),gaussian,0,0.001);MEN=double(MF)+noise;subplot(2,2,2),imshow(MEN,);title(b) 運(yùn)動(dòng)模糊+高斯噪聲后的圖像,FontSize,10)NSR=sum(noise(:).2)/sum(MEN(:).2);subplot(2,2,3),imshow(deconvwnr(MEN,PSF),); %matlab自帶維納濾波函數(shù)，若噪聲功率譜設(shè)為0，則為逆

16、濾波title(c) 逆濾波后的圖像,FontSize,10)subplot(2,2,4),imshow(deconvwnr(MEN,PSF,NSR),);%維納濾波title(d) 維納濾波后的圖像,FontSize,10)Matlab約束濾波函數(shù)：deconvreg 奇異值（奇異值（SVDSVD）偽逆濾波）偽逆濾波 ( SVD=Singular Value Decomposition)HUVu vTTiiii把退化矩陣把退化矩陣H進(jìn)行進(jìn)行SVD分解，并得到分解，并得到H的廣義逆矩陣。從而可的廣義逆矩陣。從而可以把逆濾波表示成迭代形式，通過(guò)控制迭代次數(shù)，可以在恢復(fù)以把逆濾波表示成迭代形式，通

17、過(guò)控制迭代次數(shù)，可以在恢復(fù)效果較好時(shí)停止迭代，提供了一種人機(jī)交互機(jī)會(huì)。效果較好時(shí)停止迭代，提供了一種人機(jī)交互機(jī)會(huì)?；舅枷牖舅枷肫渲衅渲蠻的列向量是的列向量是AAT的特征值，特征值， V的列向量是的列向量是ATA的特征值；的特征值； 是對(duì)是對(duì)角矩陣，對(duì)角線上的值角矩陣，對(duì)角線上的值 i稱為稱為H的奇異值。的奇異值。退化矩陣退化矩陣H的的SVD分解為分解為于是恢復(fù)圖像為于是恢復(fù)圖像為111HU Vu vTTiiii則退化矩陣則退化矩陣H的廣義逆矩陣為的廣義逆矩陣為1-1fH gu v gTiiii寫成迭代形式寫成迭代形式1(1)( )ffu v gTttttt7 7 圖象灰度校正和幾何畸變校正

18、圖象灰度校正和幾何畸變校正一、灰度校正一、灰度校正 假設(shè)理想圖象為假設(shè)理想圖象為f（x,y），由于灰度失真因），由于灰度失真因子子D（x,y）的影響，實(shí)際得到的圖象為）的影響，實(shí)際得到的圖象為g（x,y），）， g（x,y）=D（x,y）f（x,y） 灰度校正是要從畸變的圖象灰度校正是要從畸變的圖象g（x,y）中復(fù)原原）中復(fù)原原 始圖象始圖象f（x,y）。一種最直接的方法是用光密度計(jì)）。一種最直接的方法是用光密度計(jì) 測(cè)量出被拍攝景物中的某一部分區(qū)域測(cè)量出被拍攝景物中的某一部分區(qū)域S內(nèi)真實(shí)的灰內(nèi)真實(shí)的灰 （亮）度數(shù)值（亮）度數(shù)值C，而對(duì)應(yīng)的圖象灰度為，而對(duì)應(yīng)的圖象灰度為gc（x,y），）， 則則

19、 代入圖象復(fù)原方程則有代入圖象復(fù)原方程則有Sy)(x,Cy)(x,gy)D(x,c),(),(y)f(x,yxyxgCcg 三、幾何畸變校正三、幾何畸變校正 在圖象獲取或顯示過(guò)程中可能產(chǎn)生圖象的幾何失真，在圖象獲取或顯示過(guò)程中可能產(chǎn)生圖象的幾何失真，如下圖所示如下圖所示 原始網(wǎng)格圖象枕形失真圖象桶形失真圖象從上圖可以看出，幾何畸變是將無(wú)失真坐標(biāo)系中從上圖可以看出，幾何畸變是將無(wú)失真坐標(biāo)系中函函 數(shù)數(shù)f（x,y）變換到另外一個(gè)坐標(biāo)上，例如，原）變換到另外一個(gè)坐標(biāo)上，例如，原先在（先在（x,y）點(diǎn)上的象素（灰度）變化到（）點(diǎn)上的象素（灰度）變化到（u,v），），在圖象上在圖象上 反映有些位置被擠壓

20、，而另一些位置被反映有些位置被擠壓，而另一些位置被擴(kuò)張。我們希望擴(kuò)張。我們希望 找到這兩個(gè)坐標(biāo)系之間的關(guān)系。找到這兩個(gè)坐標(biāo)系之間的關(guān)系。幾何基準(zhǔn)圖像的坐標(biāo)系統(tǒng)用幾何基準(zhǔn)圖像的坐標(biāo)系統(tǒng)用(x, y)來(lái)表示來(lái)表示需要校正的圖像的坐標(biāo)系統(tǒng)用需要校正的圖像的坐標(biāo)系統(tǒng)用(x, y)表示表示設(shè)兩個(gè)圖像坐標(biāo)系統(tǒng)之間的關(guān)系用解析式表示設(shè)兩個(gè)圖像坐標(biāo)系統(tǒng)之間的關(guān)系用解析式表示10101010NiNjjiijNiNjjiijyxbyyxax通常通常h1(x,y)和和h2(x,y)用多項(xiàng)式來(lái)表示：用多項(xiàng)式來(lái)表示：),(),(21yxhyyxhx通常用線性畸變來(lái)近似較小的幾何畸變通常用線性畸變來(lái)近似較小的幾何畸變yb

21、xbbyyaxaax210210更精確一些可以用二次型來(lái)近似更精確一些可以用二次型來(lái)近似2542321025423210ybxybxbybxbbyyaxyaxayaxaax 若基準(zhǔn)圖像為若基準(zhǔn)圖像為f(x,y)，畸變圖像為，畸變圖像為g(x,y)，對(duì)于景物，對(duì)于景物上的同一個(gè)點(diǎn)，假定其灰度不變，則上的同一個(gè)點(diǎn)，假定其灰度不變，則) , (),(yxgyxfv 若已知若已知h1(x,y)和和h2(x,y)通常用線性畸變來(lái)近似較小的幾何畸變通常用線性畸變來(lái)近似較小的幾何畸變),(),(),(002001yxhyxh 零級(jí)內(nèi)插零級(jí)內(nèi)插點(diǎn)最近的),(是離),(),(),(000000畸變變圖像上的網(wǎng)y

22、xyxgyxf 雙線性插值雙線性插值(,)( x , y ),( x, y ),( x , y),( x, y) 000000001111落落在在網(wǎng)網(wǎng)格格 中 中則：則：000000000000) 1, 1( ) 1,( ) 1 (), 1)(1 ( ),() 1)(1 (),(yxyxyxyxyxgyxfv h1(x,y)和和h2(x,y)未知未知通常用已知的多對(duì)對(duì)應(yīng)點(diǎn)來(lái)確定系數(shù)通常用已知的多對(duì)對(duì)應(yīng)點(diǎn)來(lái)確定系數(shù)a, b 線性畸變線性畸變可由基準(zhǔn)圖找出三個(gè)點(diǎn)可由基準(zhǔn)圖找出三個(gè)點(diǎn)(r1, s1), (r2, s2), (r3, s3)與畸變與畸變圖像上三個(gè)點(diǎn)圖像上三個(gè)點(diǎn)(u1, v1), (u2

23、, v2), (u3, v3)一一對(duì)應(yīng)。一一對(duì)應(yīng)。ybxbbyyaxaax210210將對(duì)應(yīng)點(diǎn)代入，有：將對(duì)應(yīng)點(diǎn)代入，有：210332211321210332211321111111bbbsrsrsrvvvaaasrsrsruuu解聯(lián)立方程組，得出解聯(lián)立方程組，得出6個(gè)系數(shù)。個(gè)系數(shù)。此時(shí)此時(shí) 二次型畸變二次型畸變2542321025423210ybxybxbybxbbyyaxyaxayaxaax有有12個(gè)未知量，需要個(gè)未知量，需要6對(duì)已知對(duì)應(yīng)點(diǎn)對(duì)已知對(duì)應(yīng)點(diǎn)),(),(),(),(),(),(),(),(),(),(),(),(665544332211665544332211yxyxyxyxyx

24、yxyxyxyxyxyxyx代入上式代入上式,520266626662222222221112111621111aaayyxxyxyyxxyxyyxxyxxxxAax 同樣有同樣有Aby 解方程組，得到解方程組，得到ai, bi. 12個(gè)系數(shù)。個(gè)系數(shù)。記做向量矩陣形式記做向量矩陣形式: 通常實(shí)際應(yīng)用中，會(huì)取多余的對(duì)應(yīng)點(diǎn)對(duì)，這時(shí)通常實(shí)際應(yīng)用中，會(huì)取多余的對(duì)應(yīng)點(diǎn)對(duì)，這時(shí)A不是方陣而是高矩陣，這時(shí)矩陣的逆用廣義逆矩陣不是方陣而是高矩陣，這時(shí)矩陣的逆用廣義逆矩陣來(lái)求解。高矩陣的廣義逆矩陣為來(lái)求解。高矩陣的廣義逆矩陣為TTHHHH1)(在廣告制作和計(jì)算機(jī)動(dòng)畫中常常要使物體變形。在廣告制作和計(jì)算機(jī)動(dòng)畫中常

25、常要使物體變形。 內(nèi)插方法除了零級(jí)內(nèi)插和雙線性內(nèi)插外，還有內(nèi)插方法除了零級(jí)內(nèi)插和雙線性內(nèi)插外，還有B樣條插值和樣條插值和sinc插值函數(shù)內(nèi)插插值函數(shù)內(nèi)插wwwcsin)(sin 幾何畸變復(fù)原的一套方法也可以用于使圖像畸變幾何畸變復(fù)原的一套方法也可以用于使圖像畸變的工作中。的工作中。 Lkniinjjkikijkxvuax1200 Lkniinjjkikijkyvuby1200 最小。clearf=checkerboard(24);subplot(2,2,1),imshow(f)title(a) 原始圖像,FontSize,10)s=0.7;theta=pi/6;T=s*cos(theta) s

26、*sin(theta) 0; -s*sin(theta) s*cos(theta) 0; 0 0 1;tform=maketform(affine,T);g1=imtransform(f,tform,nearest);subplot(2,2,2),imshow(g1)title(b) 最鄰近插值變換,FontSize,10)g2=imtransform(f,tform);subplot(2,2,3),imshow(g2)title(c) 雙線性插值變換,FontSize,10)g3=imtransform(f,tform,Fillvalue,0.5);subplot(2,2,4),imshow(g3)title(e) 修改c的背景為灰色,FontSize,10)clearf=imread(Fig0515(a)(base-with-control-points).tif);subplot(2,2,1),imshow(f)title(a) 原始圖像,FontSize,10)inputpoints=83 81;450 56;43 293;249 392;436 442;outputpoints=68 66;375 47;42 286;275 434;523 532;tform=cp2tf