《數(shù)圖》第5章 圖像復(fù)原

數(shù)字圖像處理基礎(chǔ)DigitalImageProcessing第五章圖像復(fù)原1DigitalImageProcessing

圖像復(fù)原（恢復(fù)imagerestoration）－－圖像處理中經(jīng)典問(wèn)題之一，難解問(wèn)題之一。圖像復(fù)原目的－－恢復(fù)降質(zhì)圖像的本來(lái)面目。圖像復(fù)原方法－－弄清降質(zhì)原因，建立數(shù)學(xué)模型，逆降質(zhì)的過(guò)程恢復(fù)圖像。圖像降質(zhì)因素：光學(xué)系統(tǒng)的像差、衍射、畸變、攝像感光元件的非線性、成像過(guò)程的相對(duì)運(yùn)動(dòng)、大氣的湍流效應(yīng)、環(huán)境隨機(jī)噪聲等。圖像復(fù)原方法的分類(lèi)：處理域：空域、頻域；約束條件：無(wú)約束條件、有約束條件；復(fù)原系統(tǒng)：線性復(fù)原、非線性復(fù)原；…………2DigitalImageProcessing圖像復(fù)原和圖像增強(qiáng)的相似之處：改善圖像質(zhì)量。圖像增強(qiáng)－－提高圖像的“可懂度”，有較好的觀賞效果，具一定的主觀成分，沒(méi)有明確的客觀標(biāo)準(zhǔn)；圖像復(fù)原－－提高圖像的“保真度”，恢復(fù)原圖像，有更多的客觀成分，具有明確的客觀標(biāo)準(zhǔn)。3DigitalImageProcessing第1節(jié)圖像的降質(zhì)模型1.圖像降質(zhì)圖像復(fù)原處理的關(guān)鍵：建立圖像的降質(zhì)（退化）模型。常用線性、空間不變系統(tǒng)模型近似。

圖5.1圖像的退化及復(fù)原模型g(x,y)f(x,y)H[·]n(x,y)T[·]退化模型復(fù)原模型4DigitalImageProcessing2.圖像降質(zhì)模型（1）連續(xù)圖像的降質(zhì)模型不考慮噪聲時(shí)輸出圖像（5.1）其中（5.2）

（5.2）代入（5.1）式后（5.3）設(shè)H是線性算子（5.4）5DigitalImageProcessing其中

H[]是移不變算子

(5.5)線性和空間不變系統(tǒng)的退化圖像g(x,y)

(5.6)

降質(zhì)系統(tǒng)的輸出＝輸入圖像信號(hào)*系統(tǒng)沖激響應(yīng)當(dāng)沖激響應(yīng)函數(shù)h(x,y)已知時(shí)，從f(x,y)得到g(x,y)容易，從g(x,y)恢復(fù)f(x,y)并非易事。6DigitalImageProcessing考慮加性噪聲n(x,y)：

(5.7)對(duì)應(yīng)的頻率域表達(dá)式：

(5.8)

一般頻域圖像復(fù)原步驟：已知g(x,y)、h(x,y)求f(x,y)

轉(zhuǎn)化已知G(u,v)、H(u,v)求F(u,v)；g(x,y)、h(x,y)求f(x,y)

FT頻域傳遞函數(shù)H(u,v)；由G(u,v)=H(u,v)F(u,v)IFT得到復(fù)原圖像f(x,y)。7DigitalImageProcessing（2）離散圖像的降質(zhì)模型連續(xù)圖像復(fù)原推廣離散圖像復(fù)原。

空域卷積表示（離散線性空間不變系統(tǒng)）：

降質(zhì)圖像g(i,j)輸入圖像f(i,j)降質(zhì)系統(tǒng)h(i,j,m,n)加性噪聲n(i,j)

(5.9)

空域矩陣表表示：

g、f、n分別為g、f、n的堆疊列矢量

(5.10)

頻域DFT表示：

(5.11)8DigitalImageProcessing

3.圖像降質(zhì)實(shí)例（1）孔徑衍射造成的圖像降質(zhì)物平面上的點(diǎn)光源（二維沖激函數(shù)）

像平面上的光斑（系統(tǒng)沖激響應(yīng)）物平面上場(chǎng)景＝眾多點(diǎn)光源的集合

像平面上圖像＝眾多光斑的集合。點(diǎn)擴(kuò)展函數(shù)（PSF）SOd1d2點(diǎn)光源像光斑物平面像平面(x,y)圖5.2簡(jiǎn)單的透鏡成像系統(tǒng)成像透鏡(ξ,η)9在實(shí)際光學(xué)系統(tǒng)中：透鏡成像系統(tǒng)＝孔徑（aperture）成像系統(tǒng)，孔徑衍射效應(yīng)是造成圖像模糊的主要原因?？讖剑喝〈哥R的一張不透光的平面，在光軸中心開(kāi)有透光的（圓）小孔，孔中的透光率服從某種分布。衍射效應(yīng)：輸入點(diǎn)光源發(fā)出的球形發(fā)散入射光波，響應(yīng)不再是球形會(huì)聚的出射光波的單點(diǎn)圖像。圖5.3孔徑衍射造成的模糊(a)理想的PSF(b)衍射引起的PSF10DigitalImageProcessing

光瞳函數(shù)：透鏡等效孔徑所在平面上透光率的空間分布。

(5.12)光學(xué)成像的惠更斯-菲涅爾原理：對(duì)于相干光，點(diǎn)擴(kuò)展函數(shù)在幅值上就是光瞳函數(shù)的二維傅立葉變換。即：

(5.13)

進(jìn)行變量替換

(5.14)11DigitalImageProcessing相干光：降質(zhì)系統(tǒng)傳遞函數(shù)H(u,v)是h(x,y)的傅立葉變換,點(diǎn)擴(kuò)展函數(shù)h(x,y)是光瞳函數(shù)p(ξ,η)的傅立葉變換,系統(tǒng)的傳遞函數(shù)是光瞳函數(shù)的兩次傅立葉變換。

對(duì)一個(gè)函數(shù)作兩次光學(xué)傅立葉變換相當(dāng)于將此函數(shù)繞原點(diǎn)作反轉(zhuǎn)，相干光學(xué)成像系統(tǒng)的傳遞函數(shù)：

(5.15)非相干光:

各點(diǎn)光源的相位是隨機(jī)獨(dú)立的，不具有統(tǒng)一的幅度關(guān)系，只具有統(tǒng)計(jì)意義上的強(qiáng)度關(guān)系：降質(zhì)系統(tǒng)的點(diǎn)擴(kuò)展函數(shù)是相干光的點(diǎn)擴(kuò)展函數(shù)模的平方（光瞳函數(shù)FT平方）。

非相干光系統(tǒng)的點(diǎn)擴(kuò)展函數(shù)是光瞳函數(shù)的功率譜。h(x,y)H(u,v)p(ξ,η)FTFT12DigitalImageProcessing

（2）相對(duì)運(yùn)動(dòng)造成的圖像模糊景物和相機(jī)之間的相對(duì)運(yùn)動(dòng)，造成圖像的模糊，為了簡(jiǎn)單起見(jiàn)，只研究由勻速直線運(yùn)動(dòng)所造成的模糊圖像的恢復(fù)，假設(shè)：相機(jī)曝光時(shí)感光單元所產(chǎn)生的圖像退化只受相對(duì)運(yùn)動(dòng)影響之外。(a)原始圖像(b)運(yùn)動(dòng)造成的模糊圖像(c)復(fù)原后的圖像圖5.4相對(duì)運(yùn)動(dòng)造成的圖像模糊及其復(fù)原13DigitalImageProcessing設(shè)：物體f(x,y)在一平面內(nèi)運(yùn)動(dòng)，x0(t)是物體在x方向的位移，y0(t)是物體在y方向的位移，t表示運(yùn)動(dòng)的時(shí)間；感光單元的總曝光量是在快門(mén)打開(kāi)到關(guān)閉這段曝光時(shí)間T內(nèi)的積分。曝光成像后的降質(zhì)圖像為：（5.18）（5.20）（5.19）勻速直線運(yùn)動(dòng)所造成圖像模糊系統(tǒng)的傳遞函數(shù)。FT14DigitalImageProcessing再簡(jiǎn)化：物體只在x方向上作勻速直線運(yùn)動(dòng)（y0(t)=0）。圖像模糊后任意一點(diǎn)的值（5.18）式可簡(jiǎn)化為（5.21）設(shè)物體在T時(shí)間內(nèi)移動(dòng)距離為s，速度為s/T，水平方向移動(dòng)的距離為降質(zhì)系統(tǒng)的傳遞函數(shù)為： （5.22）對(duì)應(yīng)的系統(tǒng)點(diǎn)擴(kuò)展函數(shù)為：（5.23）15DigitalImageProcessing（3）大氣湍流造成的圖像降質(zhì)受大氣湍流的影響：航空?qǐng)D片、衛(wèi)星圖片、天文圖，…………簡(jiǎn)單情況：在長(zhǎng)時(shí)間作用的情況下，大氣湍流降質(zhì)圖像的系統(tǒng)傳遞函數(shù)為： （5.24）C為與湍流性質(zhì)有關(guān)的常數(shù)。16DigitalImageProcessing4.降質(zhì)模型的矩陣表示先考慮一維情況，然后再推廣到二維。（1）一維模型的矩陣表示

如f(x)為信號(hào)，h(x)為降質(zhì)系統(tǒng)的沖激響應(yīng)，則降質(zhì)系統(tǒng)的輸出g(x)為：

(5.25)不考慮噪聲上式的矩陣形式：(5.26)

式中g(shù)、f都是M維列向量，H是M×M階矩陣，矩陣H中的每一行元素均相同，每行以循環(huán)方式右移一位。

矩陣H是循環(huán)矩陣。g=Hf17DigitalImageProcessing（2）二維模型的矩陣表示圖像為f(x,y)，降質(zhì)系統(tǒng)響應(yīng)為h(x,y)，輸出降質(zhì)圖像g(x,y)為：

(5.27)用矩陣表示：(5.28)

式中g(shù)、f是MN×1維列向量，降質(zhì)矩陣H可以表示為M×M的分塊循環(huán)矩陣，其中任意一子矩陣Hj

是由降質(zhì)函數(shù)h(x,y)的第j行構(gòu)成的，也是循環(huán)矩陣。如果考慮到加性噪聲：

(5.29)用矩陣表示：(5.30)g=Hf

g=Hf＋n18DigitalImageProcessing對(duì)于實(shí)際應(yīng)用，要想從式（5.30）得出f(x,y)，其計(jì)算工作是十分困難的：

（5.30）解方程計(jì)算量十分驚人

利用循環(huán)矩陣、塊循環(huán)矩陣的性質(zhì)來(lái)對(duì)矩陣計(jì)算進(jìn)行大幅度的削減，例如將循環(huán)矩陣對(duì)角化來(lái)簡(jiǎn)化矩陣求逆運(yùn)算。

（5.30）往往是一種病態(tài)方程

對(duì)噪聲或誤差非常敏感，微小的數(shù)據(jù)誤差可以引起方程解的巨大波動(dòng)。對(duì)方程解施加一定的限制或約束，將解方程問(wèn)題轉(zhuǎn)化為求函數(shù)優(yōu)化或函數(shù)極值問(wèn)題。19DigitalImageProcessing

5.循環(huán)矩陣的對(duì)角化

利用循環(huán)矩陣的相似性定理，將降質(zhì)模型中的大矩陣進(jìn)行對(duì)角化。循環(huán)矩陣具有正交、復(fù)指數(shù)特征向量，將矩陣方程轉(zhuǎn)化成DFT的形式。龐大方程組（空間域）運(yùn)算簡(jiǎn)化頻譜運(yùn)算（頻率域）（1）循環(huán)矩陣的對(duì)角化若H、B都是M×M矩陣，若有可逆矩陣U，使得B和H互為相似矩陣H=UBU-1(5.31)

相似矩陣的特征多項(xiàng)式相同，對(duì)應(yīng)的特征值也相同。如果M階矩陣H與對(duì)角陣D相似，充要條件是H有M個(gè)線性無(wú)關(guān)的特征向量。相似矩陣D和H的特征值相同，對(duì)角陣D的對(duì)角元素就是自己的特征值。20DigitalImageProcessingM階循環(huán)矩陣H，存在M個(gè)線性無(wú)關(guān)的特征向量，可和一個(gè)對(duì)角矩陣相似。M階循環(huán)矩陣H，其特征值為一標(biāo)量集合{λ(k),k=0,1,…,M-1}，其中

(5.32)對(duì)應(yīng)的特征向量為一列矢量集合{w(k),k=0,1,…,M-1}，其中

(5.33)將上述特征值和特征向量直接代入下式驗(yàn)證，等式成立，

Hw(k)=λ(k)w(k)k=0,1,…,M-1

(5.34)

將H的M個(gè)特征向量w(k)組成一個(gè)M×M的矩陣W：

W=[w(0)w(2)…w(M-1)(5.35)只要是M×M的循環(huán)矩陣，它的特征向量組成的矩陣W都是相同的。21DigitalImageProcessing由于循環(huán)矩陣H具有M個(gè)線性無(wú)關(guān)的特征向量｛w(k)｝，可以保證由這些列向量組成的矩陣W的逆矩陣存在，即WW-1=I。參照(5.31)式的H=UBU-1

，可把H對(duì)角化：

H=WDW-1 (5.36)

D是對(duì)角陣，其對(duì)角元素d(k,k)是H的特征值λ(k)：

(5.37)

22DigitalImageProcessingH矩陣的逆矩陣為：

H-1=(WDW-1)-1=(W-1D-1W)T=WD-1W-1

(5.38）D是對(duì)角陣，D-1也是對(duì)角陣，它的每個(gè)元素為D中對(duì)應(yīng)元素的倒數(shù);W矩陣的值和H矩陣的元素?zé)o關(guān)，因此H-1比較容易獲得。H的特征值是h(x)的DFT：

h(x)形成的循環(huán)矩陣H的第k個(gè)特征值＝h(x)的DFT的第k個(gè)值H(k)。

H的相似對(duì)角矩陣D的對(duì)角元素＝序列h(x)的DFT。23DigitalImageProcessing（2）分塊循環(huán)矩陣的對(duì)角化

分塊循環(huán)矩陣H（MN×MN，M×M塊，子塊為N×N）可以對(duì)角化；其MN個(gè)特征向量組成的矩陣W，是一個(gè)的循環(huán)分塊矩陣。W的第i行第m列個(gè)子塊可表示式為：

i,m=0,1,2,…,M-1(5.39)其子塊wN為一個(gè)N×N子矩陣，其第k行第n列位置的元素為：

k,n=0,1,2,…,N-1(5.40)24DigitalImageProcessing逆矩陣W-1也是一個(gè)分塊循環(huán)矩陣。即第i行第m列個(gè)子塊表示式為：

i,m=0,1,2,…,M-1(5.41)子塊wN-1為一個(gè)N×N矩陣，其第k行n列位置的元素為：

k,n=0,1,2,…,N-1(5.42)將（5.39）至（5.42）式代入下式：WW-1

＝I，確互為逆矩陣。實(shí)際上：W是由分塊循環(huán)矩陣H的特征向量所構(gòu)成的MN×MN矩陣，其MN個(gè)特征向量是線性無(wú)關(guān)的，必定存在逆矩陣W-1。25DigitalImageProcessing分塊循環(huán)矩陣H可以寫(xiě)成：

H=WDW-1 (5.43)或 D=W-1HW(5.44)其中：矩陣D是一個(gè)MN×MN維的對(duì)角陣，由M×M個(gè)對(duì)角子矩陣組成的分塊循環(huán)矩陣，

對(duì)角子矩陣由N×N個(gè)對(duì)角元素d(k,k)組成，且為H矩陣的特征值。還可以證明，H的轉(zhuǎn)置矩陣HT可用D的復(fù)共軛D*來(lái)表示：HT=WD*W-1(5.45)26DigitalImageProcessing（3）對(duì)角化在降質(zhì)模型中的應(yīng)用1）一維情況將H=WDW-1代入（5.26）：g=Hf=WDW-1f(5.46)用W-1左乘上式兩邊：W-1g=W-1WDW-1f=DW-1f

(5.47)其中W-1是一個(gè)M×M維的矩陣，即：

(5.48)27DigitalImageProcessingg是M維列向量，W-1g也是M維列向量，其第k項(xiàng)為G(k)：即g(x)的傅立葉變換

k=0,1,2,…,M-1(5.49)同理，W-1f的第k項(xiàng)記為F(k)：

即f(x)的傅立葉變換

k=0,1,2,…,M-1(5.50)可見(jiàn)：用W-1乘g，產(chǎn)生的新向量W-1g，它的元素是g中各元素的傅立葉變換。（5.46）的g=Hf=WDW-1f中的D的主對(duì)角線元素是H的特征值λ(k)，根據(jù)（5.35）有：

k=0,1,2,…,M-1(5.51)其中H(k)是h(x)的傅立葉變換。28DigitalImageProcessing綜上，將式（5.47）簡(jiǎn)化成傅立葉變換的乘積：

(5.52)

等價(jià)的頻率域處理：(5.53)小結(jié)：

求解M個(gè)方程的方程組問(wèn)題計(jì)算三個(gè)序列的傅立葉變換問(wèn)題。W-1g=W-1WDW-1f=DW-1f29DigitalImageProcessing2）二維情況考慮噪聲項(xiàng)，g=Hf+n，用W-1左乘此式，

W-1g=W-1(Hf+n)=W-1Hf+W-1n=W-1WDW-1f+W-1n

(5.54)

化簡(jiǎn)：W-1g=DW-1f+W-1n (5.55)類(lèi)似一維情況，(5.55)式的左邊MN×1維列矢量W-1g可用下式表示：

(5.56)其中G(u,v)：

(5.57)同理，W-1f：

W-1n：u=0,1,2,…,M-1v=0,1,2,…,N-130DigitalImageProcessing類(lèi)似一維，對(duì)角陣D的任一元素為H矩陣的特征值＝降質(zhì)函數(shù)h(x,y)的DFT。 （5.58）將對(duì)角陣D的對(duì)角線元素“拉直”排列后成為：

對(duì)應(yīng)g=Hf+n，可得到

u=0,1,2,…,M-1，v=0,1,2,…,N-1(5.59)G(u,v)/F(u,v)/N(u,v)/H(u,v)分別為g(x,y)/f(x,y)/n(x,y)/h(x,y)的2D-DFT。降質(zhì)模型的龐大方程組計(jì)算簡(jiǎn)化M×N的2D-DFT計(jì)算。式中的下標(biāo)uv并不代表D矩陣的第u行第v列，而表示在D矩陣對(duì)角線上、和右乘列矢量相對(duì)應(yīng)的第（u,v）個(gè)元素。31DigitalImageProcessing第2節(jié)無(wú)約束圖像復(fù)原循環(huán)矩陣的對(duì)角化解決矩陣H求逆問(wèn)題，對(duì)噪聲的影響未予考慮，沒(méi)有對(duì)問(wèn)題的解施加一定的約束－－無(wú)約束圖像復(fù)原：

如逆濾波，運(yùn)動(dòng)模糊的消除，最小二乘方法等。1.逆濾波圖像復(fù)原

不考慮圖像噪聲，降質(zhì)圖像、降質(zhì)系統(tǒng)、原始圖像的關(guān)系為：G(u,v)=H(u,v)F(u,v)

F(u,v)=G(u,v)/H(u,v)f(x,y)=IDFT[F(u,v)]

(5.60)

H(u,v)≠0

(5.61)(5.62)復(fù)原函數(shù)1/H(u,v)和降質(zhì)函數(shù)H(u,v)作用相反:逆濾波(inversefiltering)復(fù)原。(unconstraintimagerestoration)32DigitalImageProcessing考慮噪聲G(u,v)=H(u,v)F(u,v)+N(u,v)時(shí)，復(fù)原圖像為：

(5.63)如H(u,v)出現(xiàn)零點(diǎn)，就會(huì)導(dǎo)致不定解或不穩(wěn)定解（病態(tài)性質(zhì)）：

沒(méi)有噪聲，一般也不可能精確地復(fù)原f(x,y)；考慮噪聲，噪聲項(xiàng)將被放大，使得噪聲成分在圖像復(fù)原的結(jié)果中起主導(dǎo)作用。原圖像(b)模糊圖像(c)模糊加噪圖像(d)模糊圖像逆濾波(e)模糊加噪圖像的逆濾波圖5.5逆濾波圖像復(fù)原33DigitalImageProcessing為了克服H(u,v)零點(diǎn)的影響，限制復(fù)原在頻譜坐標(biāo)離原點(diǎn)不太遠(yuǎn)的中低頻、信噪比高有限區(qū)域內(nèi)運(yùn)行；在H(u,v)零點(diǎn)的那些頻譜點(diǎn)及其附近，人為地設(shè)置H-1(u,v)的值。例如取圖像復(fù)原函數(shù)M(u,v)為：

(5.64)其中d為常數(shù)，在以d為半徑的頻域區(qū)內(nèi)，H(u,v)無(wú)零點(diǎn)，按照1/H(u,v)進(jìn)行逆濾波，在此區(qū)域以外，則不加改變。34DigitalImageProcessing

2.運(yùn)動(dòng)模糊的消除實(shí)際中：拍攝快速運(yùn)動(dòng)物體，鏡頭在曝光瞬間的偏移、運(yùn)動(dòng)的目標(biāo)……用逆濾波復(fù)原法去除由勻速直線運(yùn)動(dòng)引起的圖像模糊。（1）逆濾波復(fù)原運(yùn)動(dòng)模糊的降質(zhì)函數(shù)（5.22）為： （5.65）

對(duì)H(u,v)=0的處理：在進(jìn)行逆濾波時(shí)避免零點(diǎn)，在范圍內(nèi)進(jìn)行。35DigitalImageProcessing（2）遞推法復(fù)原在只有x方向上的運(yùn)動(dòng)（不考慮y），降質(zhì)圖像為

0≤x≤L

（5.66）要求解原圖像，需解此積分方程。經(jīng)過(guò)變量代換、上式對(duì)x求導(dǎo)數(shù)，得到：

0≤x≤L

（5.68）為了簡(jiǎn)便，設(shè)L=Ks，K為正整數(shù)，s表示在曝光期間T內(nèi)圖像中景物移動(dòng)的總距離。將圖像劃為寬度為s的K個(gè)豎條，在0～L范圍內(nèi)任意一點(diǎn)的x坐標(biāo)：

x=ms+z（5.69）

式中m為x/s的整數(shù)部分，其值為0,1,2,…,K-1，余數(shù)用z表示。sz(ms+z)msyx運(yùn)動(dòng)物體圖像f(x,y)圖5.6水平勻速運(yùn)動(dòng)降質(zhì)圖像的復(fù)原示意圖36DigitalImageProcessing將（5.69）代入（5.68）得： （5.70）

它表示原圖像的值之間是一種疊代關(guān)系。

……疊代過(guò)程的描述省略……得出只有x方向勻速直線運(yùn)動(dòng)模糊圖像復(fù)原后的表達(dá)式： （5.78）

按照上述思路，將y代入，即可得到最終復(fù)原圖像：（5.79）第m條中一點(diǎn)的原圖像值降質(zhì)圖像的導(dǎo)數(shù)值鄰近第(m-1)條中相同位置的原圖像值A(chǔ)為常數(shù)，當(dāng)K很大時(shí)接近f(x)的平均值。37DigitalImageProcessing3.無(wú)約束最小二乘方復(fù)原

圖像復(fù)原－－給定降質(zhì)圖像g，根據(jù)降質(zhì)系統(tǒng)H和噪聲n的特性，估計(jì)原圖像，使其盡量接近原圖像，兩者的誤差盡可能小。－－典型的誤差函數(shù)最小二乘方優(yōu)化問(wèn)題。降質(zhì)模型g=Hf+n，不考慮噪聲，g=Hf根據(jù)上述準(zhǔn)則，即尋找一個(gè)，使g與之偏差在最小二乘方意義上最小，即使得的范數(shù)最小，用準(zhǔn)則函數(shù)表示：

(5.80)(5.81)38DigitalImageProcessing無(wú)約束復(fù)原：在準(zhǔn)則函數(shù)求極小值的過(guò)程中，不受任何其它條件的約束。對(duì)求導(dǎo)并令其等于零：

(5.84)解得無(wú)約束最小二乘方復(fù)原圖像：

(5.85)式中(HTH)-1HT

是H矩陣的廣義逆矩陣。如H為非奇異方陣，是線性方程組g=Hf的直接解。39DigitalImageProcessing按(5.38)式的H-1=WD-1W-1，用對(duì)角化矩陣D來(lái)描述，復(fù)原圖像為：

(5.86)根據(jù)對(duì)角陣與傅立葉變換的關(guān)系：上式等價(jià)于：

(5.87)可見(jiàn)，在H為非奇異方陣、不考慮噪聲的情況下：最小二乘方準(zhǔn)則下尋找出的最優(yōu)估計(jì)圖像由g和降質(zhì)系統(tǒng)H的對(duì)角化矩陣D得出復(fù)原圖像兩者相同40DigitalImageProcessing第3節(jié)有約束圖像復(fù)原回顧：無(wú)約束圖像復(fù)原，不受任何其它條件的約束，具有相當(dāng)大的局限性。改進(jìn)：增加約束條件－－有約束圖像復(fù)原。圖像復(fù)原線性方程組是一個(gè)病態(tài)問(wèn)題，解波動(dòng)性很大。為了獲得良態(tài)解（更好的圖像復(fù)原效果），根據(jù)實(shí)際圖像復(fù)原的具體要求，利用降質(zhì)模型、噪聲、圖像的先驗(yàn)特性，對(duì)問(wèn)題的解附加一定的限制和約束，增加解的可靠性和實(shí)用性，獲得更加接近原圖像的復(fù)原效果。41DigitalImageProcessing1.有約束最小二乘方復(fù)原

有約束圖像復(fù)原：已降質(zhì)系統(tǒng)、噪聲、圖像的先驗(yàn)特性作為約束條件。如：有約束最小二乘方復(fù)原、功率譜均衡復(fù)原、平滑約束濾波復(fù)原等。

有約束最小二乘方復(fù)原：采用拉格朗日乘子法，尋找，使準(zhǔn)則函數(shù)最小：

令Q為的線性算子，要設(shè)法尋找一個(gè)最優(yōu)估計(jì)，使形式為、服從約束條件的函數(shù)最小化。求這類(lèi)問(wèn)題的最小化，用Q算子對(duì)解進(jìn)行約束，使得它的范數(shù)平方最小。拉格朗日乘子法的約束條件拉格朗日乘數(shù)此約束的物理含義：估計(jì)誤差等于噪聲。對(duì)隨機(jī)變量，“等于”是在范數(shù)平方意義上的相等。(5.88)最小二乘方準(zhǔn)則函數(shù)。42DigitalImageProcessing對(duì)求偏導(dǎo)并使其為零后解得： （5.89）式中γ=1/α，選擇不同的Q，可得到不同的有約束復(fù)原方法。當(dāng)Q=I，對(duì)解不產(chǎn)生約束，即對(duì)解不進(jìn)行線性運(yùn)算，解得： （5.90）

當(dāng)γ=0時(shí)，退化為無(wú)約束情況的圖像復(fù)原，

前（5.85）

可見(jiàn)無(wú)約束最小二乘方復(fù)原是有約束最小二乘方復(fù)原的特殊情況。在空間域上求解(5.89)或(5.90)，由于維數(shù)巨大，很難直接求解，通常轉(zhuǎn)換到頻率域上進(jìn)行處理。43DigitalImageProcessing

2.維納濾波復(fù)

維納濾波復(fù)原：約束算子Q=Rf-1/2Rn1/2

Rf

=E{ffT}（堆疊圖像f的相關(guān)矩陣）Rn=E{nnT

}（堆疊噪聲n的相關(guān)矩陣）

（5.91）

Rf的第i

行第j

列元素是E{fi

fj}，表示圖像f的第i個(gè)和第j個(gè)元素的相關(guān)，Rn的第i

行第j

列元素是E{ni

nj}，表示噪聲n的第i個(gè)和第j個(gè)元素的相關(guān)。Rf

和Rn

都是實(shí)對(duì)稱矩陣。

在圖像中，鄰近像素點(diǎn)高度相關(guān)，距離較遠(yuǎn)的像素其相關(guān)性較弱。典型的相關(guān)矩陣在主對(duì)角線方向呈一條非零元素帶，其它區(qū)域多為零值或接近零值。44DigitalImageProcessing默認(rèn)圖像和噪聲的平穩(wěn)性，可近似將Rf和Rn

用塊循環(huán)矩陣表示：Rf

=WDfW-1 Rn

=WDnW-1

（5.92）根據(jù)塊循環(huán)矩陣對(duì)角化的原理，列矢量表示受約束后的：（5.94）……得到復(fù)原后的圖像： （5.100）對(duì)應(yīng)的頻域表達(dá)式：（5.101）45DigitalImageProcessing下面討論式（5.101）的幾種情況：①（5.101）式：有信噪比最大化約束條件下的最小二乘方復(fù)原濾波器。②如γ=1，（5.101）式為維納濾波器，如γ為變數(shù)，則為參變維納濾波器。③無(wú)噪聲時(shí)，Sn(u,v)=0。式（5.101）退化成逆濾波器。④當(dāng)無(wú)信號(hào)時(shí)，這時(shí)濾波器都完全截止，阻止噪聲通過(guò)。⑤如果無(wú)噪聲的統(tǒng)計(jì)性質(zhì)，一般可假設(shè)為白噪聲＝常數(shù)功率譜密度。⑥當(dāng)Sf(u,v)和Sn(u,v)統(tǒng)計(jì)性質(zhì)未知時(shí)，式（5.101）可以用下式近似： （5.103）46DigitalImageProcessing圖5.7維納濾波圖像復(fù)原一例，復(fù)原后的峰值信噪比（PSNR）提高4dB以上。(a)原始圖像(b)加噪模糊圖像(22dB)(c)復(fù)原后圖像(26dB)圖5.7維納濾波圖像復(fù)原47DigitalImageProcessing3.功率譜均衡復(fù)原

（powerspectrumequalization）復(fù)原原則：估計(jì)圖像的功率譜＝原圖像的功率譜相等 （5.104）由此導(dǎo)出功率譜均衡復(fù)原濾波器：

（5.105） 根據(jù)圖像退化的降質(zhì)模型和功率譜的定義可知：（5.106）估計(jì)圖像功率譜原圖像功率譜復(fù)原函數(shù)降質(zhì)圖像功率譜48DigitalImageProcessing將(5.106)式代入(5.105)式，可得: （5.107） （5.108）（5.109）(5.106)(5.105)功率譜均衡vs維納濾波相似之處：需要已知信號(hào)和噪聲的功率譜；無(wú)噪聲時(shí)，簡(jiǎn)化為逆濾波；無(wú)信號(hào)時(shí)，完全截止。不同之處：功率譜均衡濾波器無(wú)零點(diǎn)問(wèn)題。49DigitalImageProcessing4.平滑約束復(fù)原(smoothnessconstraints)自然場(chǎng)景的圖像數(shù)據(jù)間的變化是比較平緩的，突然的變化較少，因此，要求所估計(jì)的圖像數(shù)據(jù)是平滑過(guò)渡的，是施加了平滑約束的。

平滑約束最小二乘方復(fù)原：使得估計(jì)圖像的二階導(dǎo)數(shù)（平滑的標(biāo)志）最小。只需要知道有關(guān)噪聲的均值和方差等先驗(yàn)知識(shí)。約束算子矩陣Q：選擇拉普拉斯算子，避免一般導(dǎo)數(shù)的方向性。差分式近似：

(5.110)等效卷積模板：(5.111)

50DigitalImageProcessing準(zhǔn)則函數(shù)中Q是MN×MN維的分塊循環(huán)矩陣（平滑矩陣）：（5.112）為了滿足矩陣運(yùn)算要求，將p(x,y)矩陣補(bǔ)零后形成N×N維pe(x,y)矩陣，上式中每個(gè)子矩陣qj

是pe(x,y)矩陣的第j行組成的N×N維的循環(huán)矩陣，即（5.113)51DigitalImageProcessing重寫(xiě)（5.89）式最優(yōu)解的結(jié)果：（5.114）按照對(duì)角化方法：（5.115）上式兩邊左乘以W-1，得到（5.116）復(fù)原圖像的頻譜表示：（5.117）

形式上與維納濾波器有些相似，但這里只需對(duì)噪聲均值和方差的知識(shí)；在平滑約束中，獲得的最優(yōu)解是經(jīng)拉普拉斯濾波后其能量之和最小的一幅圖像；拉普拉斯濾波的結(jié)果表示圖像的跳變分量（高頻分量）大??；選擇跳變分量最少，符合實(shí)際圖像具有高度相關(guān)、以中低頻能量為主的特性。52DigitalImageProcessing第4節(jié)非線性圖像復(fù)原實(shí)際中：圖像降質(zhì)系統(tǒng)并非線性，噪聲并非加性，線性復(fù)原只是一種近似。近年來(lái)，非線性圖像復(fù)原的方法獲得了良好的效果。1.最大后驗(yàn)概率復(fù)原

(MAP，MaximumAPosterior)將原始圖像f(x,y)和被觀測(cè)到的降質(zhì)圖像g(x,y)都看成二維隨機(jī)場(chǎng)，在已知g(x,y)的前提下，對(duì)f(x,y)進(jìn)行種種估計(jì):求出后驗(yàn)概率密度函數(shù)P(f(x,y)|g(x,y))，當(dāng)P(f(x,y)|g(x,y))為最大值時(shí)所對(duì)應(yīng)的就代表在已知降質(zhì)圖像g(x,y)時(shí)最有可能的原始圖像f(x,y)，也就是是f(x,y)的最大后驗(yàn)估計(jì)。 根據(jù)貝葉斯準(zhǔn)則求P(f|g)的最大值：

(5.118)降質(zhì)圖像的分布P(g)獨(dú)立于最優(yōu)化問(wèn)題，可用常數(shù)C表示。53DigitalImageProcessing忽略常數(shù)，對(duì)式（5.118）取對(duì)數(shù)后再求偏導(dǎo)使其為零，解出f即滿足要求：

(5.119)考慮先前的降質(zhì)模型：g=Hf+n (5.120)

設(shè)：H是服從N維高斯分布的平穩(wěn)隨機(jī)場(chǎng)，只需均值和方差參數(shù)就可以確定。噪聲n滿足零均值的N維高斯分布：

(5.121)f的分布已知，g的隨機(jī)性是由n決定的，P(g|f)和P(n)具有相同分布：

(5.122)上式對(duì)f求偏導(dǎo)可得：(5.123)n的協(xié)方差矩陣54DigitalImageProcessing設(shè)圖像像素f的概率密度函數(shù)為正態(tài)分布：

(5.124)對(duì)f求偏導(dǎo)：

(5.125)把(5.123)、(5.125)代入(5.119)：

(5.126)如事先通過(guò)使f的均值為0，上式可解得：

(5.127)這是關(guān)于的線性方程，可用疊代或其它方法求解得到恢復(fù)圖像。(5.119)(5.123)55DigitalImageProcessing

如果噪聲n為N維獨(dú)立高斯分布隨機(jī)場(chǎng)，(5.127)式中的協(xié)方差矩陣為單位矩陣，將不起作用，維納濾波中的約束算子相當(dāng)于這里的，在這種情況下MAP方法和最小二乘方復(fù)原方法是一致的。注意：如(5.118)式中得不到f的分布，或簡(jiǎn)單地將f作為均勻分布，則此時(shí)MAP復(fù)原就成為最大似然估計(jì)(ML，MaximumLikehood)復(fù)原：

(5.128)(5.127)(5.89)MAP結(jié)果比較

有約束最小二乘方結(jié)果56DigitalImageProcessing

2.最大熵復(fù)原最大熵復(fù)原：是在圖像復(fù)原過(guò)程中加最大熵約束，要求在所有可行解中，圖像熵最大那一個(gè)作為最終的圖像解。

優(yōu)點(diǎn)：不需要對(duì)圖像先驗(yàn)知識(shí)做更多假設(shè)，在抑制噪聲和恢復(fù)圖像細(xì)節(jié)之間取得較好的平衡。信息論：熵的大小是表征隨機(jī)變量集合（信源）的隨機(jī)程度的統(tǒng)計(jì)平均值；信源為等概率分布式，其熵值最大。如果將N×N大小的非負(fù)圖像f(i,j)像素值看作是一個(gè)隨機(jī)變量集合，任一像素的灰度值表示這一點(diǎn)隨機(jī)變量的概率，仿效隨機(jī)變量熵，定義圖像熵Hf：

(5.129)這里的圖像熵和以往的圖像信源熵是不同的。57DigitalImageProcessing信源熵：隨機(jī)集合的熵值越大，說(shuō)明各個(gè)隨機(jī)變量的概率大小越趨同，差別越小。圖像熵：圖像的熵值越大，說(shuō)明各個(gè)像素的灰度值越趨同，差別越小，圖像越平滑。本質(zhì)上：圖像最大熵的約束就是圖像平滑的約束。噪聲熵Hn：保證圖像的非負(fù)性，定義等效噪聲n’(i,j)=n(i,j)+B(5.130)

(5.131)圖像最大熵復(fù)原問(wèn)題等價(jià)求圖像熵和噪聲熵加權(quán)熵的極大值問(wèn)題（Frieden方法）。加權(quán)熵H：

(5.132)為加權(quán)系數(shù)，用以調(diào)整對(duì)噪聲熵的強(qiáng)調(diào)，值大，噪聲熵Hn對(duì)總量H的貢獻(xiàn)也大，對(duì)噪聲的平滑要求加強(qiáng)，結(jié)果對(duì)噪聲的平滑作用要遠(yuǎn)遠(yuǎn)超過(guò)對(duì)圖像的平滑，以免對(duì)圖像過(guò)分平滑。最大的噪聲負(fù)值58DigitalImageProcessing在用拉各朗日法對(duì)（5.132）式求極值的過(guò)程中還要增加具體的約束條件。第一個(gè)約束：圖像的灰度值之和為一常數(shù)

(5.133)第二個(gè)約束：降質(zhì)圖像g(p,q)為h(p,q)和原