基于MATLAB的圖像恢復(fù)算法研究

1、PAGE PAGE 28中北大學(xué)課 程 設(shè) 計 說 明 書學(xué)生姓名： 學(xué) 號： 學(xué)生姓名： 學(xué) 號： 學(xué)生姓名： 學(xué) 號： 學(xué) 院： 信息商務(wù)學(xué)院 專 業(yè)： 電子信息工程 題 目： 信息處理綜合實踐: 基于MATLAB的圖像恢復(fù)算法研究 指導(dǎo)教師： 職稱: 年 月 日中北大學(xué)課程設(shè)計任務(wù)書 13/14 學(xué)年第 一 學(xué)期學(xué) 院： 信息商務(wù)學(xué)院 專 業(yè)： 電子信息工程 學(xué) 生 姓 名： 學(xué) 號： 學(xué) 生 姓 名： 學(xué) 號： 學(xué) 生 姓 名： 學(xué) 號: 課程設(shè)計題目： 信息處理綜合實踐: 于MATLAB的圖像恢復(fù)算法研究 起 迄 日 期： 課程設(shè)計地點： 電子信息科學(xué)與技術(shù)專業(yè)實驗室 指 導(dǎo) 教 師

2、： 系 主 任： 下達任務(wù)書日期: 年月日課 程 設(shè) 計 任 務(wù) 書1設(shè)計目的：1、通過本課程設(shè)計的學(xué)習(xí)，學(xué)生將復(fù)習(xí)所學(xué)的專業(yè)知識，使課堂學(xué)習(xí)的理論知識應(yīng)用于實踐，通過本課程設(shè)計的實踐使學(xué)生具有一定的實踐操作能力；2、掌握Matlab使用方法，能熟練運用該軟件設(shè)計并完成相應(yīng)的信息處理；3、通過圖像處理實踐的課程設(shè)計，掌握設(shè)計圖像處理軟件系統(tǒng)的思維方法和基本開發(fā)過程。2設(shè)計內(nèi)容和要求（包括原始數(shù)據(jù)、技術(shù)參數(shù)、條件、設(shè)計要求等）：(1)創(chuàng)建一個仿真運動/均值模糊PSF來模糊一幅圖像（圖像自選）；(2)針對退化設(shè)計出復(fù)原濾波器，對退化圖像進行復(fù)原(復(fù)原的方法自定)；(3)對退化圖像進行復(fù)原，顯示復(fù)原

3、前后圖像，對復(fù)原結(jié)果進行分析，并評價復(fù)原算法；(4)要求每位學(xué)生進行查閱相關(guān)資料，并寫出自己的報告。注意每個學(xué)生的報告要有所側(cè)重，寫出自己所做的內(nèi)容。3設(shè)計工作任務(wù)及工作量的要求包括課程設(shè)計計算說明書(論文)、圖紙、實物樣品等：每個同學(xué)獨立完成自己的任務(wù)，每人寫一份設(shè)計報告，在課程設(shè)計論文中寫明自己設(shè)計的部分，給出設(shè)計結(jié)果。課 程 設(shè) 計 任 務(wù) 書4主要參考文獻：阮秋琦等.數(shù)字圖像處理（第三版）.北京：電子工業(yè)出版社.20112. 岡薩雷斯等.數(shù)字圖像處理（MATLAB版）.北京：電子工業(yè)出版社.20013. 楊杰等.數(shù)字圖像處理及MATLAB實現(xiàn)：學(xué)習(xí)與實驗指導(dǎo).北京：電子工業(yè)出版社.20

4、104. 劉衛(wèi)國等.MATLAB程序設(shè)計與應(yīng)用.北京：高等教育出版社.20065. 許國根等.模式識別與智能計算的MATLAB實現(xiàn).北京：北京航空航天大學(xué)出版社.20125設(shè)計成果形式及要求：畢業(yè)設(shè)計說明書仿真結(jié)果6工作計劃及進度：2013年12月16日 12月19日：查資料；12月19日 12月24日：在指導(dǎo)教師指導(dǎo)下設(shè)計方案；12月25日 12月27日：撰寫課程設(shè)計說明書； 12月27日：答辯系主任審查意見： 簽字： 年 月 日目 錄 TOC o 1-3 h z u HYPERLINK l _Toc375851342 摘要： PAGEREF _Toc375851342 h 5 HYPERL

5、INK l _Toc375851343 1圖像復(fù)原的概念 PAGEREF _Toc375851343 h 6 HYPERLINK l _Toc375851344 1.1圖像復(fù)原的定義 PAGEREF _Toc375851344 h 6 HYPERLINK l _Toc375851345 1.2 圖象恢復(fù)與圖象增強的異同 PAGEREF _Toc375851345 h 6 HYPERLINK l _Toc375851346 1.3 圖象退化的原因 PAGEREF _Toc375851346 h 6 HYPERLINK l _Toc375851347 1.4 維納濾波的研究歷史 PAGEREF _

6、Toc375851347 h 6 HYPERLINK l _Toc375851348 1.5圖象退化舉例 PAGEREF _Toc375851348 h 7 HYPERLINK l _Toc375851349 2退化模型 PAGEREF _Toc375851349 h 7 HYPERLINK l _Toc375851350 2.1圖象退化模型概述 PAGEREF _Toc375851350 h 7 HYPERLINK l _Toc375851351 2.2連續(xù)函數(shù)退化模型 PAGEREF _Toc375851351 h 8 HYPERLINK l _Toc375851352 2.3離散函數(shù)退化

7、模型 PAGEREF _Toc375851352 h 8 HYPERLINK l _Toc375851353 3.圖象復(fù)原技術(shù) PAGEREF _Toc375851353 h 9 HYPERLINK l _Toc375851354 3.1無約束恢復(fù) PAGEREF _Toc375851354 h 9 HYPERLINK l _Toc375851355 3.2 逆濾波 PAGEREF _Toc375851355 h 9 HYPERLINK l _Toc375851356 3.3 維納（Wiener）濾波器基本原理 PAGEREF _Toc375851356 h 9 HYPERLINK l _To

8、c375851357 3.4維納濾波復(fù)原法 PAGEREF _Toc375851357 h 11 HYPERLINK l _Toc375851358 3.5圖像復(fù)原例圖 PAGEREF _Toc375851358 h 12 HYPERLINK l _Toc375851359 4圖像復(fù)原的MATLAB實現(xiàn)實例 PAGEREF _Toc375851359 h 13 HYPERLINK l _Toc375851360 5.結(jié)束語 PAGEREF _Toc375851360 h 13 HYPERLINK l _Toc375851361 參考文獻： PAGEREF _Toc375851361 h 14

9、HYPERLINK l _Toc375851362 附錄： PAGEREF _Toc375851362 h 14 HYPERLINK l _Toc375851363 （1）.維納濾波復(fù)原源代碼： PAGEREF _Toc375851363 h 14 HYPERLINK l _Toc375851364 （2）.規(guī)則化濾波復(fù)原程序源代碼： PAGEREF _Toc375851364 h 14 HYPERLINK l _Toc375851365 （3）.Lucy-Richardson復(fù)原濾波源代碼： PAGEREF _Toc375851365 h 15 HYPERLINK l _Toc3758513

10、66 （4）.盲目去卷積復(fù)原源代碼： PAGEREF _Toc375851366 h 15摘要：圖像復(fù)原是圖象處理的一個重要課題。圖像復(fù)原也稱圖象恢復(fù)，是圖象處理中的一大類技術(shù)。它的主要目的是改善給定的圖像質(zhì)量。當(dāng)給定了一幅退化了的或者受到噪聲污染了的圖像后，利用退化現(xiàn)象的某種先驗知識來重建或恢復(fù)原有圖像是復(fù)原處理的基本過程?？赡艿耐嘶泄鈱W(xué)系統(tǒng)中的衍射，傳感器非線性畸變，光學(xué)系統(tǒng)的像差，攝影膠片的非線性，大氣湍流的擾動效應(yīng)，圖像運動造成的模糊及幾何畸變等等。噪聲干擾可以由電子成像系統(tǒng)傳感器、信號傳輸 過程或者膠片顆粒性造成。各種退化圖像的復(fù)原都可歸結(jié)為一種過程，具體地說就是把退化模型化，并且

11、采用相反的過程進行處理，以便恢復(fù)出原圖像。文章介紹了圖象退化的原因，幾種常用的圖像濾波復(fù)原技術(shù)，以及用MATLAB實現(xiàn)圖像復(fù)原的方法。1圖像復(fù)原的概念1.1圖像復(fù)原的定義圖像復(fù)原也稱圖象恢復(fù)，是圖象處理中的一大類技術(shù)。所謂圖像復(fù)原，是指去除或減輕在獲取數(shù)字圖像過程中發(fā)生的圖像質(zhì)量下降（退化）這些退化包括由光學(xué)系統(tǒng)、運動等等造成圖像的模糊，以及源自電路和光度學(xué)因素的噪聲。圖像復(fù)原的目標(biāo)是對退化的圖像進行處理，使它趨向于復(fù)原成沒有退化的理想圖像。成像過程的每一個環(huán)節(jié)（透鏡，感光片，數(shù)字化等等）都會引起退化。在進行圖像復(fù)原時，既可以用連續(xù)數(shù)學(xué)，也可以用離散數(shù)學(xué)進行處理。其次，處理既可在空間域，也可在

12、頻域進行。1.2 圖象恢復(fù)與圖象增強的異同相同點：改進輸入圖像的視覺質(zhì)量 。不同點：圖象增強目的是取得較好的視覺結(jié)果(不考慮退化原因)； 圖象恢復(fù)根據(jù)相應(yīng)的退化模型和知識重建或恢復(fù)原始的圖像(考慮退化原因)。1.3 圖象退化的原因 圖象退化指由場景得到的圖像沒能完全地反映場景的真實內(nèi)容，產(chǎn)生了失真等問題。其原因是多方面的。如： 透鏡象差/色差 聚焦不準（失焦，限制了圖像銳度） 模糊（限制頻譜寬度） 噪聲（是一個統(tǒng)計過程）抖動（機械、電子）1.4 維納濾波的研究歷史維納是著名的數(shù)學(xué)家，后來被譽為信息理論家。維納的著作不僅是一個很好的創(chuàng)見，而且具有結(jié)合工程的實際意義，是線性濾波理論研究的一個重要的

13、開端. 在第二次世界大戰(zhàn)中，由于雷達的發(fā)明以及防空炮火控制的任務(wù)，把大量有修養(yǎng)的數(shù)學(xué)家和物理學(xué)家都動員到信息科學(xué)這個研究領(lǐng)域中來了，這個時候人們活躍于這個領(lǐng)域，并有許多重大的科學(xué)創(chuàng)造。數(shù)學(xué)家維納對于濾波理論的研究成果，就是這時候重大的科學(xué)創(chuàng)見之一。 通訊與控制中的濾波問題，指的是從獲得的信號與干擾中盡可能地濾除干擾，分離出所期望的信號，或者說，是通過對一系列帶有誤差的實際測量數(shù)據(jù)的處理，得出所期望數(shù)據(jù)的估計值川。維納的工作是從研究處在統(tǒng)計平衡的時間序列開始的，維納證明:在一定條件下，處在統(tǒng)計平衡的時間序列的時間平均等于相平均。維納正是基于這點提出了他著名的濾波和預(yù)測理論。濾波問題就是盡可能地恢

14、復(fù)一個被噪聲干擾了的信號的問題。實質(zhì)上，就是預(yù)測一個被噪聲干擾了的時間序列的問題，因此，濾波問題也可以視為一個預(yù)測問題。數(shù)學(xué)上講，預(yù)測就是從一個時間序列的過去的數(shù)據(jù)估算整個序列的統(tǒng)計參數(shù)。 工程上的濾波問題也是理論上的一類統(tǒng)計估計問題，最佳線形濾波是最佳線性估計的方法之一，在最佳估計中最小均方誤差估計是最有現(xiàn)實意義的。估計理論的課題是眾多的，最小均方誤差估計只是估計理論的一個小的分支。然而，它卻是最重要又最富有實際意義的一個分支，對系統(tǒng)所加的線性條件起初是為了簡化理論分析，非線性濾波問題是在理論處理上比線性濾波問題要困難和復(fù)雜的多，但是后來證明:在一定條件下，在最小均方誤差準則下得到的最佳線性

15、系統(tǒng)是所有系統(tǒng)中的最佳者。 近代濾波理論的發(fā)展對于信息科學(xué)的發(fā)展是有重大貢獻的，它概括了通訊與控制中信息過濾的統(tǒng)計本質(zhì)。這是由于濾波理論與通訊和控制中的許多課題有密切的聯(lián)系，從而賦予了濾波理論以極大的生命力，濾波理論本來是一個小的研究領(lǐng)域，但是它聯(lián)系著許多大的廣泛的研究領(lǐng)域，因此它的價值己經(jīng)超出了它起源時自身的價值，也就是它能夠繼續(xù)活躍地向前發(fā)展的保證。 幾十年來濾波理論已經(jīng)發(fā)展成了一個廣闊的研究領(lǐng)域，可以有許多不同的方法來介紹它的內(nèi)容，有的可以選擇不同的重點。本文主要是關(guān)于維納濾波的，介紹維納濾波的基本概念以及講其維納濾波的應(yīng)用。 從數(shù)學(xué)的觀點來說濾波理論是統(tǒng)計學(xué)中的估計理論的一個重要分支，

16、從工程的觀點來看它又是系統(tǒng)工程研究的一個重要組成部分。1.5圖象退化舉例如圖1所示是兩個圖象退化的例子。圖1 退化圖像與原始圖像2退化模型2.1圖象退化模型概述圖像復(fù)原處理的關(guān)鍵問題在于建立退化模型。在用數(shù)學(xué)方法描述圖像時，它的最普遍的數(shù)學(xué)表達式為這樣一個表達式可以代表一幅活動的、彩色的立體圖像。當(dāng)研究的是靜止的、單色的、平面的圖像時，則其數(shù)學(xué)表達式就簡化為基于這樣的數(shù)學(xué)表達式，可建立如圖2所示的退化模型。由圖2的模型可見，一幅純凈的圖像是由于通過了一個系統(tǒng)H及加性噪聲而使其退化為一幅圖像的。圖2 圖像退化模型圖像復(fù)原可以看成是一個估計過程。如果已經(jīng)給出了退化圖像并估計出系統(tǒng)參數(shù)H，從而可近似

17、地恢復(fù)。這里，是一種統(tǒng)計性質(zhì)的噪聲信息。當(dāng)然，為了對處理結(jié)果做出某種最佳的估計，一般應(yīng)首先明確一個質(zhì)量標(biāo)準。根據(jù)圖像的退化模型及復(fù)原的基本過程可見，復(fù)原處理的關(guān)鍵在于對系統(tǒng)H的基本了解。就一般而言，系統(tǒng)是某些元件或部件以某種方式構(gòu)造而成的整體。退化模型可分為連續(xù)函數(shù)退化模型和離散函數(shù)退化模型。2.2連續(xù)函數(shù)退化模型假定系統(tǒng)H對坐標(biāo)為(,)處的沖激函數(shù)(x-,y-)的沖激響應(yīng)為h(x,y,)，則此式說明，如果系統(tǒng)H對沖激函數(shù)的響應(yīng)為已知，則對任意輸入的響應(yīng)可用上式求得，即，線性系統(tǒng)H完全可以由沖激響應(yīng)來表征。圖像中沖激響應(yīng)也稱為點擴散函數(shù)。在有噪音的情況下：2.3離散函數(shù)退化模型對和進行均勻取樣

18、后，就可引伸出離散函數(shù)的退化模型。用一維的來說明。如果f (x)和h(x)周期分別A和B的序列，為避免卷積周期重疊需要對它們進行周期擴展為周期為M A + B 1。 f(x) 0 x A-1 h(x) 0 x B-1fe(x)= he(x)=0 A-1 x M-1 0 B-10.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,

19、Facecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(u

20、int8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth低通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift

21、(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth高通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r

22、(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯低通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=

23、freqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);fi

24、guresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯高通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=

25、Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形低通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqr

26、t(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong)

27、; 梯形高通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplo

28、t(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 用其他方法編寫的理想低通、理想高通、Butterworth低通、同態(tài)濾波程序理想低通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=flo

29、or(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=1; else h=0; end y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 理想高通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g)