第12章圖像復原

第12章圖像復原12.1圖像退化與復原12.2非約束復原12.3最小二乘類約束復原12.4非線性復原方法12.5其他圖像復原技術12.6編程實例12.1圖像退化與復原

12.1.1圖像降質(zhì)的數(shù)學模型

圖像復原處理的關鍵問題在于建立退化模型。輸入圖像f(x，y)經(jīng)過某個退化系統(tǒng)后的輸出是一幅退化的圖像。為了討論方便，把噪聲引起的退化即噪聲對圖像的影響一般作為加性噪聲考慮，這也與許多實際應用情況一致。如圖像數(shù)字化時的量化噪聲、隨機噪聲等就可以作為加性噪聲，即使不是加性噪聲而是乘性噪聲，也可以用對數(shù)方式轉(zhuǎn)化為相加形式。

原始圖像f(x，y)經(jīng)過一個退化算子或退化系統(tǒng)H(x，y)的作用，并且和噪聲n(x，y)進行疊加，形成退化后的圖像g(x，y)。圖12-1表示退化過程輸入和輸出的關系。圖中H(x，y)概括了退化系統(tǒng)的物理過程，就是所要尋找的退化數(shù)學模型。圖12-1圖像的退化模型數(shù)字圖像的圖像恢復問題可看做：根據(jù)退化圖像g(x，y)和退化算子H(x，y)的形式，沿著反向過程去求解原始圖像f(x，y)，或者說是逆向?qū)ふ以紙D像的最佳近似估計。圖像退化的模型過程可以用數(shù)學表達式寫成如下形式：(12-1)在這里n(x，y)是一種統(tǒng)計性質(zhì)的信息。在實際應用中，往往假設噪聲是白噪聲，即它的頻譜密度為常數(shù)，并且與圖像不相關。

在圖像復原處理中，盡管非線性、時變和空間變化的系統(tǒng)模型更具有普遍性和準確性，更與復雜的退化環(huán)境相接近，但它給實際處理工作帶來巨大的困難，常常找不到解或者很難用計算機來處理。因此，在圖像復原處理中，往往用線性系統(tǒng)和空間不變系統(tǒng)模型來加以近似。這種近似的優(yōu)點使得線性系統(tǒng)中的許多理論可直接用于解決圖像復原問題，同時不失可用性。下面介紹連續(xù)圖像退化的數(shù)學模型。

一幅連續(xù)圖像f(x，y)可以看做是由一系列點源組成的。因此，f(x，y)可以通過點源函數(shù)的卷積來表示。即(12-2)式中：δ為點源函數(shù)，表示空間上的點脈沖。在不考慮噪聲的一般情況下，連續(xù)圖像經(jīng)過退化系統(tǒng)H后的輸出為(12-3)把式(12-2)代入式(12-3)得(12-4)在線性和空間不變系統(tǒng)的情況下，退化算子H具有如下性質(zhì)。性質(zhì)1線性。設f1(x，y)和f2(x，y)為兩幅輸入圖像，k1和k2為常數(shù)，則(12-5)由性質(zhì)1還可推出下面兩個結(jié)論：(1)當k1=k2=1時，式(12-5)變?yōu)?12-6)(2)如f2(x，y)=0，則(12-7)性質(zhì)2空間不變性。如果對任意f(x，y)以及a和b，有(12-8)對于線性空間不變系統(tǒng)，輸入圖像經(jīng)退化后的輸出為(12-9)式中，h(x-α，y-β)為該退化系統(tǒng)的點擴展函數(shù)，或叫系統(tǒng)的沖激響應函數(shù)。它表示系統(tǒng)對坐標為(α，β)處的沖激函數(shù)δ(x-α，y-β)的響應。也就是說，只要系統(tǒng)對沖激函數(shù)的響應為已知，那么就可以清楚圖像退化是如何形成的。因為對于任一輸入f(α，β)的響應，都可以通過上式計算出來。此時，退化系統(tǒng)的輸出就是輸入圖像信號f(x，y)與點擴展函數(shù)h(x，y)的卷積:(12-10)圖像退化除了受到成像系統(tǒng)本身的影響外，有時還要受到噪聲的影響。假設噪聲n(x，y)是加性白噪聲，這時上式可寫成：(12-11)在頻域上，式(12-11)可以寫成：(12-12)上式中，G(u，v)、F(u，v)、N(u，v)分別是退化圖像g(x，y)、原圖像f(x，y)、噪聲信號n(x，y)的傅立葉變換；H(u，v)是系統(tǒng)的點沖激響應函數(shù)h(x，y)的傅立葉變換，稱為系統(tǒng)在頻率域上的傳遞函數(shù)。式(12-11)和式(12-12)就是連續(xù)函數(shù)的退化模型?？梢?，圖像復原實際上就是已知g(x，y)求f(x，y)的問題或已知G(u，v)求F(u，v)的問題，它們的不同之處在于一個是在空間域，一個是在頻域。

顯然，進行圖像復原的關鍵問題是尋找降質(zhì)系統(tǒng)在空間域上的沖激響應函數(shù)h(x，y)，或者降質(zhì)系統(tǒng)在頻率域上的傳遞函數(shù)H(u，v)。一般來說，傳遞函數(shù)比較容易求得。因此，在進行圖像復原之前，應設法求得完全的或近似的降質(zhì)系統(tǒng)傳遞函數(shù)。要想得到h(x，y)，只需對H(u，v)求傅立葉逆變換即可。12.1.2離散圖像退化的數(shù)學模型

1.一維離散退化模型

設f(x)為具有A個采樣值的離散輸入函數(shù)，h(x)為具有B個采樣值的退化系統(tǒng)的沖激響應函數(shù)，則經(jīng)退化系統(tǒng)后的離散輸出函數(shù)g(x)為輸入f(x)和沖激響應h(x)的卷積，即為了避免上述卷積所產(chǎn)生的各個周期重疊(設每個采樣函數(shù)的周期為M)，分別對f(x)和h(x)用添0延伸的方法擴展成周期M=A+B-1的周期函數(shù)，即(12-13)輸出為(12-14)式中：x=0，1，2，…，M-1。因為fe(x)和he(x)已擴展成周期函數(shù)，故ge(x)也是周期性函數(shù)，用矩陣表示為(12-15)因為he(x)的周期為M，所以he(x)=he(x+M)，即

M×M階矩陣H可寫為(12-16)上式寫成更簡潔的形式：(12-17)式中，g、f都是M維列向量；H是M×M階矩陣，矩陣中的每一行元素均相同，只是每行以循環(huán)方式右移一位，因此矩陣H是循環(huán)矩陣。循環(huán)矩陣相加或相乘得到的還是循環(huán)矩陣。2.二維離散模型

設輸入的數(shù)字圖像f(x，y)大小為A×B，點擴展函數(shù)h(x，y)被均勻采樣為C×D大小。為避免交疊誤差，仍用添0擴展的方法，將它們擴展成M=A+C-1和N=B+D-1個元素的周期函數(shù)：(12-18)則輸出的降質(zhì)數(shù)字圖像為(12-19)式中：x=0，1，2，…，M-1；y=0，1，2，…，N-1。式(12-19)的二維離散退化模型同樣可以用式(12-17)所示的矩陣形式表示，即

g=Hf式中：g、f是MN×1維列向量；H是MN×MN維矩陣。用矩陣形式表示二維離散退化模型的方法是將g(x，y)和f(x，y)中的元素排成列向量。(12-20)(12-21)

Hi為子矩陣，大小為N×N，即H矩陣是由M×M個大小為N×N的子矩陣組成的，稱為分塊循環(huán)矩陣。分塊矩陣是由延拓函數(shù)he(x，y)的第j行構(gòu)成的，構(gòu)成方法如下：(12-22)若把噪聲考慮進去，則離散圖像退化模型為(12-23)寫成矩陣形式為

g=Hf+n (12-24)上述線性空間不變退化模型表明，在給定g(x，y)，并且知道退化系統(tǒng)的點擴展函數(shù)h(x，y)和噪聲分布n(x，y)的情況下，可估計出原始圖像f(x，y)。

假設圖像大小為MN=512×512=262144，其相應矩陣H的元素個數(shù)也為262144個，這意味著要解出f(x，y)，需要解262144個聯(lián)立方程組，其計算量十分驚人。考慮到矩陣H為循環(huán)矩陣，因此可利用循環(huán)矩陣的性質(zhì)簡化運算，限于篇幅，本書不做討論。12.2非約束復原

非約束復原是指在已知退化圖像g的情況下，根據(jù)對退化系統(tǒng)H和n的一些了解或假設，估計出原始圖像，使得某種事先所確定的誤差準則為最小。12.2.1逆濾波

由式(12-24)可得

n=g-Hf(12-25)

逆濾波法是指在對n沒有先驗知識的情況下，可以依據(jù)這樣的最優(yōu)準則，即尋找一個，使得H在最小二乘方誤差的意義下最接近g，即要使n的?；蚍稊?shù)(norm)最?。?12-26)求式(12-26)的最小值：

如果我們在求最小值的過程中不做任何約束，則稱這種復原為非約束復原。

由極值條件：(12-27)(12-28)解出為

對上式作傅立葉變換，得

可見，如果知道g(x，y)和h(x，y)，也就知道了G(u，v)和H(u，v)。根據(jù)上式，即可得出F(u，v)，再經(jīng)過反傅立葉變換就能求出f(x，y)。

逆濾波是最早應用于數(shù)字圖像復原的一種方法，用此方法處理過由漫游者、探索者等衛(wèi)星探索發(fā)射得到的圖像。(12-29)(12-30)12.2.2非約束圖像復原的病態(tài)性質(zhì)

由式(12-30)進行圖像復原時，由于H(u，v)在分母上，若在(u，v)平面上某引起點或區(qū)域H(u，v)很小或等于零，即出現(xiàn)了零點，就會導致不穩(wěn)定解。因此，即使沒有噪聲，一般也不可能精確地復原f(x，y)。如果考慮噪聲項N(x，y)，則出現(xiàn)零點時，噪聲項將被放大，零點的影響將會更大，對復原的結(jié)果起主導地位，這就是無約束圖像復原模型的病態(tài)性質(zhì)。它意味著退化圖像中小的噪聲干擾在H(u，v)取得很小值的那些頻譜上將對恢復圖像產(chǎn)生很大的影響。由簡單的光學分析知道，在超出光學系統(tǒng)的繞射極限時，H(u，v)將很小或等于零，因此對多數(shù)圖像直接采用逆濾波復原會遇到上述求解方程的病態(tài)性。為了克服這種不穩(wěn)定性，一方面可利用后面所講的有約束圖像復原；另一方面，可利用噪聲一般在高頻范圍，其衰減速度較慢，而信號的頻譜隨頻率升高下降較快的性質(zhì)，在復原時，只限制在頻譜坐標距原點不太遠的有限區(qū)域內(nèi)運行，而且關心的也是信噪比高的那些頻率位置。Nathan在用逆濾波圖像復原時采

用的是限定恢復轉(zhuǎn)移函數(shù)最大值的方法，其傳遞函數(shù)H(u，v)和恢復轉(zhuǎn)移函數(shù)M(x，y)如圖12-2所示。圖12-2逆濾波復原實際上，為了避免H(u，v)值太小，一種改進方法是在H(u，v)=0的那些頻譜點及其附近，人為地設置H-1(u，v)的值，使得在這些頻譜點附近N(u，v)/H(u，v)不會對

產(chǎn)生太大的影響。圖12-3給出了H(u，v)、H-1(u，v)應用這種改進的濾波特性或恢復轉(zhuǎn)移函數(shù)的一維波形，從中可以看出它與正常濾波的差別。圖12-3逆濾波器零點的影響及其改進另一種改進是考慮到退化系統(tǒng)的傳遞函數(shù)H(u，v)帶寬比噪聲的帶寬要窄得多，其頻率特性具有低通性質(zhì)，因此取恢復轉(zhuǎn)移函數(shù)M(u，v)為

其中，ω0的選取原則是將H(u，v)為零的點除去。這種方法的缺點是復原后圖像的振鈴效果較明顯。(12-31)12.3最小二乘類約束復原

非約束復原是指除了使準則函數(shù)最小外，再沒有其他的約束條件。因此它只需了解降質(zhì)系統(tǒng)的傳遞函數(shù)或點擴展函數(shù)，就能利用如前所述的方法進行復原。但是由于傳遞函數(shù)存在病態(tài)問題，復原只能局限在靠近原點的有限區(qū)域內(nèi)進行，使非約束圖像復原具有相當大的局限性。最小二乘類約束復原是指除了要求了解關于退化系統(tǒng)的傳遞函數(shù)之外，還需要知道某些噪聲的統(tǒng)計特性或噪聲與圖像的某些相關情況。根據(jù)所了解的噪聲的先驗知識的不同，采用不同的約束條件，從而得到不同的圖像復原技術。在最小二乘類約束復原中，要設法尋找一個最優(yōu)估計，使得形式為的函數(shù)最小化。求這類問題的最小化，常采用拉格朗日乘子算法。也就是說，要尋找一個，使得準則函數(shù)：(12-32)為最小。上式中，Q為f[TX(]的線性算子；α為一常數(shù)，稱為拉格朗日乘子。對式(12-32)求導：求解得到：(12-33)上式中γ=1/α，這個常數(shù)必須調(diào)整到約束被滿足為止。求解式(12-33)的關鍵就是如何選用一個合適的變換矩陣Q。Q的形式不同，就可得到不同類型的有約束的最小二乘類圖像復原方法。如果用圖像f和噪聲的相關矩陣Rf和Rn表示Q。就可以得到維納濾波復原方法。如選用拉普拉斯算子形式，即使某個函數(shù)的二階導數(shù)最小，就可推導出有約束最小平方恢復方法。12.3.1維納濾波

在一般情況下，圖像信號可近似地認為是平穩(wěn)隨機過程，維納濾波將原始圖像f和對原始圖像的估計作為隨機變量。假設Rf和Rn為f和n的自相關矩陣，其定義為(12-34)式中：E{·}代表數(shù)學期望運算。

Rf和Rn均為實對稱矩陣，在大多數(shù)圖像中，鄰近的像素點是高度相關的，而距離較遠的像素其相關性卻較弱。通常，f和n的元素之間的相關不會延伸到20～30個像素的距離之外。因此，一般來說，自相關矩陣在主對角線附近有一個非零元素帶，而在右上角和左上角的區(qū)域內(nèi)將為零值。如果像素之間的相關是像素之間距離的函數(shù)，而不是它們位置的函數(shù)，可將Rf和Rn近似為分塊循環(huán)矩陣。因而，用循環(huán)矩陣的對角化，可寫成：(12-35)式中：W為一個MN×MN矩陣，包含M×M個N×N的塊。M、N含義見二維離散模型部分。

W的第(i，m)個分塊為

(12-36)其中WN為一個N×N矩陣，其第(k，n)個位置的元素為

式(12-35)中，A和B的元素分別為Rf和Rn中的自相關元素的傅立葉變換。這些自相關元素的傅立葉變換被分別定義為fe(x，y)和ne(x，y)的譜密度Sf(u，v)和Sn(u，v)。定義QTQ=R-1fRn，代入式(12-33)，得：(12-37)進一步可推導出：(12-38)式中：D是對角陣，D*為D的共軛矩陣。D的對角元素與

he(x,y)中的傅立葉變換有關：對式(12-38)再進行矩陣變換假設M=N，則：(12-39)u，v分別為0，1，2，…，N-1，式中|H(u,v)|2=H*(u,v)H(u,v)。對式(12-39)作如下分析：

(1)如果γ=1，則稱之為維納濾波器。注意，當γ=1時，并不是在約束條件下得到的最佳解，即并不一定滿足

。若γ為變數(shù)，此式為參變維納濾波器。

使用參變維納濾波器時，H(u，v)由點擴展函數(shù)確定，而當噪聲是白噪聲時，Sn(u，v)為常數(shù)，可通過計算一幅噪聲圖像的功率譜Sg(u，v)求解。由于Sg(u,v)=|H(u,v)|2

Sf(u,v)+Sn(u,v)，所以Sf(u，v)可通過本式求得。(2)當無噪聲影響時，Sn(u，v)=0，稱之為理想的逆濾波器。逆濾波器可看成是維納濾波器的一種特殊情況。

(3)如果不知道噪聲的統(tǒng)計性質(zhì)，也就是Sf(u，v)和Sn(u，v)未知，則式(12-39)可以用下式近似：式中：K表示噪聲與信號的頻譜密度之比。12.3.2約束最小平方濾波

約束最小平方復原是一種以平滑度為基礎的圖像復原方法。如前所述，在進行圖像恢復計算時，由于退化算子矩陣H［·］的病態(tài)性質(zhì)，多數(shù)在零點附近的數(shù)值起伏過大，使得復原后的圖像產(chǎn)生了多余的噪聲和邊緣。約束最小平方濾波仍然以最小二乘方濾波復原公式(12-37)為基礎，通過選擇合理的Q，并優(yōu)化‖Qf‖2，從而去掉被恢復圖像的這種尖銳部分，即增加圖像的平滑性。我們知道，圖像增強的拉普拉斯算，它具有突出邊緣的作用，而則恢復了圖像的平滑性。因此，在進行圖像恢復時可將拉普拉斯算子作為約束?，F(xiàn)在的問題是如何將其表示成‖Qf‖2的形式，以便使用式(12-37)。

在離散情況下，拉普拉斯算子可用下面的差分運算實現(xiàn)：(12-40)利用f(x，y)與下面的模板算子進行卷積可實現(xiàn)上面的運算：(12-41)在離散卷積的過程中，可利用延伸f(x，y)和p(x，y)來避免交疊誤差。延伸后的函數(shù)為pc(x，y)。建立分塊循環(huán)矩陣，將平滑準則表示為矩陣形式：(12-42)式(12-42)中每個子矩陣Cj是pe(x，y)的第j行組成的N×N循環(huán)矩陣。即Cj表示如下：(12-43)根據(jù)循環(huán)矩陣的對角化可知，利用前述的矩陣W進行對角化，即(12-44)式中：E為對角矩陣，其元素為(12-45)

E(k，i)是C中元素Pe(x,y)的二維傅立葉變換。并且，可以將寫成fTCTCf，定義Q=C，則fTCTCf=‖Qf‖2。

如果要求約束條件‖g-Hf‖=‖n‖2得到滿足，在Q=C時，有：

(12-46)式(12-46)兩邊同乘以W-1，得：(12-47)式中：D*為D的共軛矩陣。所以：(12-48)式中：u,v分別為0，1，…，N-1，而且|H(u,v)|2=H*(u,v)H(u,v)。本濾波器也稱為最小平方濾波器。12.4非線性復原方法

12.4.1最大后驗復原

最大后驗復原是一種統(tǒng)計方法，它把原圖像f(x，y)和退化圖像g(x，y)都作為隨機場，在已知g(x，y)的前提下，求出后驗條件概率密度函數(shù)p(f(x，y)/g(x，y))。若

使下式最大：(12-49)則就代表已知退化圖像g(x，y)時，最可能的原始圖像f(x，y)。這種圖像復原方法稱之為最大后驗圖像復原方法。最大后驗圖像復原法把圖像看做是非平穩(wěn)隨機場，把圖像模型表示成一個平穩(wěn)隨機過程對于一個不平穩(wěn)的均值作零均值Gauss起伏，可得出求解迭代序列：(12-50)式中：k為迭代次數(shù)；和分別為f和n的方差的倒數(shù)；是隨空間而變的均值，它是一個常數(shù)，但要經(jīng)過多次迭代才能收斂到最后的解。12.4.2最大熵復原

最大熵復原方法通過最大化某種反映圖像平滑性的準則函數(shù)來作約束條件，以解決圖像復原中反向濾波法存在的病態(tài)問題。

首先簡單介紹一下熵的概念。

熵的定義為(12-51)式中p(x)為隨機變量x的概率密度。對于離散信號，熵的定義為(12-52)熵是表征隨機變量集合的隨機程度的量度。當所有隨機變量等可能性，也就是說P1=P2=…=Pm時熵最大，為H=lnM。由于概率P(k)介于0和1之間，因此最大熵的范圍在0～lnM之間，H不可能出現(xiàn)負值。在二維數(shù)字圖像中，熵的定義為(12-53)最大熵復原的原理是將f(x，y)寫成隨機變量的統(tǒng)計模型，然后在一定的約束條件下，找出用隨機變量形式表示的熵的表達式，運用求極大值的方法，求得最優(yōu)估計解。最大熵復原的含義是對的最大平滑估計。最大熵復原常用Friend和Burg兩種方法。這兩種方法的基本原理相同，這里僅介紹Friend法。首先定義一幅大小為M×N的圖像f(x，y)，顯然f(x，y)非負。

圖像的總能量E和熵分別為(12-54)(12-55)類似地可定義噪聲的熵Hn：式中：n′(x，y)=n(x，y)+B，B為最大噪聲負值?；謴途褪窃跐M足式(12-54)和圖像退化模型的約束條件下，使恢復后的圖像熵和噪聲熵達到最大。熵通常取決于f的形狀，當圖像具有均勻的灰度時最大。因此，用最大熵恢復圖像具有某種平滑性。引入拉格朗日(Lagrange)函數(shù)：

(12-57)式中：λmn(m，n=1，2，…，N)和β是拉格朗日乘子；ρ是加權(quán)因子，表示Hf和Hn相互之間的權(quán)重。和分別表示f(x，y)和n′(x，y)的估計值，則有：(12-58)(12-59)把式(12-57)分別代入式(12-58)、式(12-59)，可得：(12-60)(12-61)并且和滿足下列約束條件(12-62)(12-63)式(12-60)為圖像恢復函數(shù)。把式(12-60)和式(12-61)代入式(12-62)和式(12-63)可得(N2+1)個方程。由此聯(lián)立方程組可解得(N2+1)個未知數(shù)解λmn(m，n=1，2，…，N)，上述方程組可求得的值。12.4.3投影復原

投影復原法是用代數(shù)方程組來描述線性和非線性退化系統(tǒng)的：

g(x,y)=D［f(x,y)］+n(x,y)

其中：f(x，y)是原始圖像；g(x，y)是退化圖像；n(x，y)是系統(tǒng)噪聲；D是退化算子，表示對圖像進行某種運算。

在使用投影復原法進行圖像復原時，引進一些先驗信息附加的約束條件，可改善圖像復原效果。12.5其他圖像復原技術

12.5.1幾何畸變校正

數(shù)字圖像在獲取過程中，由于成像系統(tǒng)的非線性，成像后的圖像與原景物圖像相比會產(chǎn)生比例失調(diào)，甚至扭曲，我們把這類圖像退化現(xiàn)象稱之為幾何畸變。典型的幾何畸變?nèi)鐖D12-4所示。

幾何畸變校正一般要對失真的圖像進行精確的幾何校正。通常是先確定一幅圖像為基準，然后去校正另一幅圖像的幾何形狀。因此，幾何畸變校正一般分兩步來做：第一步是進行圖像空間坐標的變換；第二步是重新確定在校正空間各像素點的取值。圖12-4幾種典型的幾何畸變1.空間幾何坐標變換

按照一幅標準圖像g(u，v)或一組基準點去校正另一幅幾何失真圖像f(x，y)，稱之為空間幾何坐標變換。根據(jù)兩幅圖像的一些已知對應點對(也稱為控制點對)，建立起函數(shù)關系式，將失真圖像的坐標系(x，y)變換到標準圖像坐標系(u，v)，從而實現(xiàn)失真圖像按標準圖像的幾何位置校正，使f(x，y)中的每一像素點都可在g(u，v)中找到對應像素點。常采用三角形線性法實現(xiàn)空間幾何坐標變換。2.三角形線性法

圖像的幾何失真雖然是非線性的，但在一個局部小區(qū)域內(nèi)可近似認為是線性