第五章 圖像復(fù)原與重建(2)

1、數(shù)字圖像處理數(shù)字圖像處理 主講人：杜宏博第五章第五章 圖像復(fù)原與重建圖像復(fù)原與重建5.1 圖像退化圖像退化/復(fù)原處理的模型復(fù)原處理的模型5.2 噪聲模型噪聲模型5.3 空間濾波去噪空間濾波去噪5.4 頻域?yàn)V波去噪頻域?yàn)V波去噪5.5 退化函數(shù)建模退化函數(shù)建模5.6 圖像復(fù)原的方法圖像復(fù)原的方法 直接逆濾波直接逆濾波 維納濾波維納濾波5.7 圖像投影重建圖像投影重建周期噪聲的模型是二維正弦波，通過(guò)帶阻、帶通和陷波濾波器可以被有效去除。理想帶阻濾波器的表達(dá)式理想帶阻濾波器的表達(dá)式: :00001,( , )2( , )0,( , )221,( , )2WD u vDWWH u vDD u vDWD

2、u vD5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 帶阻濾波器帶阻濾波器n n階的巴特沃思帶阻濾波器階的巴特沃思帶阻濾波器22201( , )( , )1( , )nH u vD u v WD u vD高斯帶阻濾波器高斯帶阻濾波器2220( , )12( , )( , )1Du vDD u v WH u ve 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 帶阻濾波器帶阻濾波器理想帶阻濾波器理想帶阻濾波器巴特沃思帶阻濾波器巴特沃思帶阻濾波器高斯帶阻濾波器高斯帶阻濾波器5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲帶阻濾波器帶阻濾波器(a) 被正弦噪聲污染的圖像被正弦噪聲污染的圖像 (

3、b) 圖圖(a)的頻譜的頻譜(c) 巴特沃思帶阻濾波器巴特沃思帶阻濾波器 (d) 濾波效果圖濾波效果圖5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲( , )1( , )bpbrHu vHu v 帶通濾波器帶通濾波器帶通濾波器執(zhí)行與帶阻濾波器相反的操作。帶通濾波器執(zhí)行與帶阻濾波器相反的操作。( , )( , ):bpbrHu vHu v帶通濾波器的傳遞函數(shù)可根據(jù)相應(yīng)的帶阻濾波器的傳遞函數(shù)得到不直接使用，損失大量不直接使用，損失大量圖像細(xì)節(jié)。圖像細(xì)節(jié)?？衫脦V波器提取噪可利用帶通濾波器提取噪聲模式。聲模式。5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲陷波濾波器陷波濾波器阻止阻止(或通過(guò)

4、或通過(guò))事先定義的中心頻率鄰域內(nèi)的頻率。事先定義的中心頻率鄰域內(nèi)的頻率。理想陷波濾波器理想陷波濾波器巴特沃思陷波濾波器巴特沃思陷波濾波器高斯陷波濾波器高斯陷波濾波器由于傅立葉變換由于傅立葉變換是對(duì)稱的是對(duì)稱的,因此因此陷波濾波器必須陷波濾波器必須以關(guān)于原點(diǎn)對(duì)稱以關(guān)于原點(diǎn)對(duì)稱的形式出現(xiàn)。的形式出現(xiàn)。5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲陷波濾波器陷波濾波器00000,(,)(,)Du vuv半徑為中心在且在對(duì)稱的理想陷波濾波器的傳遞函數(shù)10200( , )( , )( , )1D u vDD u vDH u v或其他22 1/210022 1/2200( , )(/2)(/2) ( ,

5、 )(/2)(/2) D u vuMuvNvD u vuMuvNv其中5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 陷波濾波器陷波濾波器2012:1( , )1( , )/,nnH u vDD u vDu v階數(shù)為 的巴特沃思陷波帶阻濾波器的傳遞函數(shù)為1220( , )/,12:( , )1D u vDu vDH u ve 高斯陷波帶阻濾波器的傳遞函數(shù)為 還可以得到另一種陷波濾波器還可以得到另一種陷波濾波器,它能通過(guò)它能通過(guò)(而不是阻止而不是阻止)包含在陷波區(qū)的頻率包含在陷波區(qū)的頻率.陷波區(qū)域的形狀可以是任意的陷波區(qū)域的形狀可以是任意的(如矩形如矩形)。5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波

6、降低周期噪聲圖像退化模型：圖像退化模型：5.5 退化函數(shù)建模退化函數(shù)建模),(),(*),(),(yxyxfyxhyxg),(),(),(),(vuNvuFvuHvuG退化系統(tǒng)一般情況下是：線性，位置不變的退化系統(tǒng)退化系統(tǒng)一般情況下是：線性，位置不變的退化系統(tǒng)（1 1）線性：）線性：yxfbHyxfaHyxbfyxafH,2121（2 2）位置不變性：對(duì)任意）位置不變性：對(duì)任意, yxf有有yxgyxfH,對(duì)于線性位置不變退化，圖像復(fù)原其實(shí)就是一個(gè)圖像反卷積對(duì)于線性位置不變退化，圖像復(fù)原其實(shí)就是一個(gè)圖像反卷積的過(guò)程的過(guò)程 圖像觀察估計(jì)法圖像觀察估計(jì)法給定一幅退化圖像，但沒(méi)有退化函數(shù) H 的知識(shí)

7、，那么估計(jì)該函數(shù)的方法之一就是收集圖像自身的信息： 尋找簡(jiǎn)單結(jié)構(gòu)的子圖像 尋找受噪聲影響小的子圖像5.5 退化函數(shù)建模退化函數(shù)建模估計(jì)退化系統(tǒng)模型的三種方法估計(jì)退化系統(tǒng)模型的三種方法構(gòu)造一個(gè)估計(jì)圖像，它與觀察的子圖像有相同大小和特性 表示觀察子圖像， 表示構(gòu)造的子圖像 和 為對(duì)應(yīng)的傅立葉變換。( , )( , )?/( , )SssH u vGu vF u v),(yxgs),(yxfs),(vuGs),(vuFs假設(shè)空間不變的，由 推導(dǎo)出完全函數(shù) ),(vuHs),(vuH5.5 退化函數(shù)建模退化函數(shù)建模 圖像試驗(yàn)估計(jì)法圖像試驗(yàn)估計(jì)法 使用與被退化圖像設(shè)備相似的裝置，并得到一個(gè)脈沖的沖激響應(yīng)

8、，可以進(jìn)行較準(zhǔn)確的退化估計(jì)：( , )( , )G u vH u vA一個(gè)脈沖點(diǎn)一個(gè)脈沖點(diǎn)成像系統(tǒng)成像系統(tǒng)H 此處A是一個(gè)沖激的傅立葉變換，表示沖擊強(qiáng)度，為一常數(shù)。 右圖為一個(gè)放大的亮脈沖以及退化的沖激。( , )g x y退化圖像退化圖像 模型估計(jì)法模型估計(jì)法 建立退化模型，考慮引起退化的環(huán)境因素。建立退化模型，考慮引起退化的環(huán)境因素。22 5/6()( , )k uvH u ve 例如：例如：Hufnagel Hufnagel 等等 Stanley Stanley 的退化模型就是基于大氣湍的退化模型就是基于大氣湍流的物理特性而提出來(lái)的，其中流的物理特性而提出來(lái)的，其中k k為常數(shù)，與湍流特

9、性相關(guān)。為常數(shù)，與湍流特性相關(guān)。( (除了指數(shù)除了指數(shù)5/65/6，該公式與高斯低通濾波形式相同，該公式與高斯低通濾波形式相同.).)5.5 退化函數(shù)建模退化函數(shù)建模 模型估計(jì)法模型估計(jì)法5.5 退化函數(shù)建模退化函數(shù)建模大氣湍流模型模擬退化模糊一幅圖像：大氣湍流模型模擬退化模糊一幅圖像：劇烈湍流劇烈湍流(k=0.0025)(k=0.0025)中等湍流中等湍流(k=0.001)(k=0.001)輕微湍流輕微湍流(k=0.00025)(k=0.00025)可忽略可忽略的湍流的湍流22 5/6()( , )k uvH u ve如果已知系統(tǒng)的傳遞函數(shù) ，則根據(jù) vuFvuHvuG,vuH,vuHvuG

10、vuF,可得復(fù)原圖像的譜，經(jīng)傅氏逆變換即可得到復(fù)原圖像在忽略噪聲的影響，退化模型的傅氏變換為實(shí)際應(yīng)用時(shí)存在病態(tài)的問(wèn)題，即在 H(u,v) 等于零或非常小的數(shù)值點(diǎn)上， 將變成無(wú)窮大或非常大的數(shù)。 ),(vuF- 這就是逆濾波復(fù)原法5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波 vuNvuFvuHvuG,vuHvuNvuFvuHvuNvuHvuGvuF,),(,系統(tǒng)中存在噪聲時(shí)退化模型的傅立葉變換為：寫成逆濾波復(fù)原的方式：1）即使知道退化函數(shù)，也不能準(zhǔn)確復(fù)原圖像，因?yàn)樵肼暫瘮?shù) N(u,v) 是一個(gè)隨機(jī)函數(shù)，其傅里葉變換未知。2）如果退化是零或非常小的值，噪聲即使數(shù)值很小，但 N(u,v)/H(

11、u,v) 之比 (上式第二項(xiàng)) 可能非常大，很容易錯(cuò)估 的值。),(vuF12 ()( , )( , )( , )( , )jux vyf x yf x yN u v Hu vedudv 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波解決退化是零或非常小的值的途徑：限制濾波的頻率，使其接近原點(diǎn)值。 在頻率平面離原點(diǎn)較遠(yuǎn)的地方，H(u,v)數(shù)值較小或?yàn)榱?，因此圖像復(fù)原在原點(diǎn)周圍的有限區(qū)域內(nèi)進(jìn)行，即將退化圖像的傅立葉譜限制在沒(méi)出現(xiàn)零點(diǎn)而且數(shù)值又不是太小的有限范圍內(nèi)，即通過(guò)將頻率限制為接近原點(diǎn)分析，減少了遇到零值的幾率。 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波劇烈湍流劇烈湍流( (k

12、 k=0.0025)=0.0025)大氣湍流模型模擬退化模糊一幅圖像大氣湍流模型模擬退化模糊一幅圖像可忽略的湍流可忽略的湍流652222/)/()/(),(NvMukevuH對(duì)退化函數(shù)對(duì)退化函數(shù)H H( (u u, ,v v) )進(jìn)行精確取反并進(jìn)行逆濾波，結(jié)果如下圖。進(jìn)行精確取反并進(jìn)行逆濾波，結(jié)果如下圖。5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波全頻直接全頻直接逆濾波復(fù)原逆濾波復(fù)原半徑為半徑為4040時(shí)截止時(shí)截止H H半徑為半徑為7070時(shí)截止時(shí)截止H H半徑為半徑為8585時(shí)截止時(shí)截止H H結(jié)果表明：噪聲明顯影響結(jié)果表明：噪聲明顯影響了圖像復(fù)原結(jié)果，一般直了圖像復(fù)原結(jié)果，一般直接逆濾

13、波效果較差。接逆濾波效果較差。劇烈湍流圖劇烈湍流圖( (k k=0.0025)=0.0025)5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波最小均方誤差復(fù)原法最小均方誤差復(fù)原法 -Wiener-Wiener濾波復(fù)原濾波復(fù)原目標(biāo)：目標(biāo)： 尋找一個(gè)濾波器，使得復(fù)原后圖像尋找一個(gè)濾波器，使得復(fù)原后圖像 與原與原始圖像始圖像 的均方誤差最小。的均方誤差最小。 逆濾波沒(méi)有清楚說(shuō)明如何處理噪聲！逆濾波沒(méi)有清楚說(shuō)明如何處理噪聲！),(yxf)(22ffEe誤差度量：誤差度量： 現(xiàn)討論一種濾波復(fù)原法現(xiàn)討論一種濾波復(fù)原法-Wiener-Wiener濾波復(fù)原：濾波復(fù)原： 綜合考慮退化函數(shù)和噪聲統(tǒng)計(jì)特征。綜合考

14、慮退化函數(shù)和噪聲統(tǒng)計(jì)特征。E期望值。期望值。),(yxf因此維納濾波復(fù)原又稱為最小均方誤差復(fù)原。因此維納濾波復(fù)原又稱為最小均方誤差復(fù)原。min),(),(2yxfyxfE5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原),(),(/ ),(),(),(),(),(vuGvuSvuSvuHvuHvuHvuFf221f 誤差函數(shù)的最小值在頻域里可以通過(guò)近似圖像 的傅里葉變換來(lái)計(jì)算：葉變換為復(fù)原近似圖像的傅里換為退化圖像的傅里葉變換為退化函數(shù)的傅里葉變其中：),(),(),(vuFvuGvuH為未退化圖像的功率譜為噪聲的功率譜22),(),(),(),(vuFvuSvuNvuSf)

15、,(),(),(*vuHvuHvuH2維納濾波器維納濾波器5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原(2) (2) 未退化圖像的功率譜難以知道，可用下式近似表示：未退化圖像的功率譜難以知道，可用下式近似表示：KvuHvuHvuHvuHw221),(),(),(),(1) (1) 如果噪聲為如果噪聲為 0 0，其功率譜消失，維納濾波就退化為逆濾波。，其功率譜消失，維納濾波就退化為逆濾波。討論：討論：式中式中 K K 是根據(jù)信噪比的某種先驗(yàn)知識(shí)確定的常數(shù)。是根據(jù)信噪比的某種先驗(yàn)知識(shí)確定的常數(shù)。),(/ ),(),(),(),(),(vuSvuSvuHvuHvuHvuHfw2

16、21維納濾波復(fù)原：維納濾波復(fù)原：維納濾波需要假定下述條件成立維納濾波需要假定下述條件成立( (或近似成立或近似成立) )： 系統(tǒng)為線性、空間不變；系統(tǒng)為線性、空間不變；(1)(1) 退化圖像、原始圖像和噪聲都是均勻隨機(jī)場(chǎng)，噪聲的均退化圖像、原始圖像和噪聲都是均勻隨機(jī)場(chǎng)，噪聲的均值為零，且與圖像不相關(guān)。值為零，且與圖像不相關(guān)。5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原 維納濾波復(fù)原與逆濾波復(fù)原的比較全頻逆濾波半徑受限逆濾波維納濾波復(fù)原 (交互選擇K) 維納濾波的缺點(diǎn)： 未退化圖像和噪聲的功率譜必須是已知的； 功率比(信噪比)常數(shù)K 的估計(jì)一般還是沒(méi)有合適的解。5.6 圖像

17、復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原維納濾波器的維納濾波器的matlab實(shí)現(xiàn)實(shí)現(xiàn)deconvwnr ：Deblur image using Wiener filter SyntaxJ = deconvwnr(I,PSF)J = deconvwnr(I,PSF,NSR)J = deconvwnr(I,PSF,NCORR,ICORR)其中：I是退化圖像 PSF系統(tǒng)函數(shù)（點(diǎn)擴(kuò)散函數(shù)） NSR信噪比 NCORR：噪聲的自相關(guān)函數(shù) ICORR：退化圖像的自相關(guān)函數(shù) J：復(fù)原圖像5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wie

18、ner濾波復(fù)原濾波復(fù)原維納濾波和逆濾波復(fù)原案例：維納濾波和逆濾波復(fù)原案例：clc; clear; close clc; clear; close allall; ;f = double( imread(f = double( imread(cameraman.tifcameraman.tif););subplot(231); imshow(f,);subplot(231); imshow(f,);title(title(orginal clean imageorginal clean image); );% generate the degrade function% generate the

19、 degrade functionPSF = fspecial(PSF = fspecial(motionmotion,7,45);,7,45);subplot(232); imshow(PSF,);subplot(232); imshow(PSF,);title(title(Point spread functionPoint spread function); ); % using the PSF to degrade image% using the PSF to degrade imagegb = imfilter(f,PSF,gb = imfilter(f,PSF,circularc

20、ircular); );subplot(233); imshow(gb,);subplot(233); imshow(gb,);title(title(Blurred image caused by motionBlurred image caused by motion); ); 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原維納濾波和逆濾波復(fù)原案例：維納濾波和逆濾波復(fù)原案例：% add noise to the degraded image% add noise to the degraded imagenoise = imnoise(zeros(size(f),noi

21、se = imnoise(zeros(size(f),gaussiangaussian,0,0.1);,0,0.1);g = gb + noise;g = gb + noise;subplot(234);imshow(g,)subplot(234);imshow(g,)title(title(Blurred image with noiseBlurred image with noise) ) % inverse filtering% inverse filteringfr1 = deconvwnr(g,PSF);fr1 = deconvwnr(g,PSF);subplot(235);imsh

22、ow(fr1,)subplot(235);imshow(fr1,)title(title(inverse filtering resultinverse filtering result) ) 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原維納濾波和逆濾波復(fù)原案例：維納濾波和逆濾波復(fù)原案例：% wiener filtering% wiener filteringSn = abs(fft2(noise).2; % noise power spectrumnA = sum(Sn(:)/prod(size(noise); % noise average powerSf = abs

23、(fft2(f).2; % image power spectrumfA = sum(Sf(:)/prod(size(f); % image average powerR = nA/fA; % signal to noise ratiofr2 = deconvwnr(g,PSF,R);fr2 = deconvwnr(g,PSF,R);subplot(236);imshow(fr2,)subplot(236);imshow(fr2,)title(title(wiener filtering resultwiener filtering result) ) 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-W

24、iener濾波復(fù)原濾波復(fù)原維納濾波和逆濾波復(fù)原案例：維納濾波和逆濾波復(fù)原案例：5.6 圖像投影重建圖像投影重建 概念：投影重建一般指利用物體的多個(gè)（軸向概念：投影重建一般指利用物體的多個(gè)（軸向）投影圖像重建目標(biāo)圖像的過(guò)程。它是一類特）投影圖像重建目標(biāo)圖像的過(guò)程。它是一類特殊的圖像處理方法，輸入的是一系列的投影圖殊的圖像處理方法，輸入的是一系列的投影圖，輸出是重建圖。，輸出是重建圖。 通過(guò)投影重建可以直接的看到原來(lái)被投影的物通過(guò)投影重建可以直接的看到原來(lái)被投影的物體的某種特性的空間分布，比直觀觀測(cè)投影圖體的某種特性的空間分布，比直觀觀測(cè)投影圖要直觀的多。要直觀的多。5.6 圖像投影重建圖像投影重

25、建 Radon變換對(duì)f(x,y)的Radon變換g(t, )定義為沿由t和定義的直線l的線積分。,( , ),( cossin)kkjg tf x ydlf x yxyt dxdy 5.6 圖像投影重建圖像投影重建 Radon變換(, ),( coscos)jjg tfx yxytdxdy Radon 變換揭示了函數(shù)和投影之間的關(guān)系，若函數(shù)為f (x, y)，則不同角度下的投影可寫為原理：原理：“斷層平面中某一點(diǎn)的密度值可看作這一平面內(nèi)所有經(jīng)過(guò)斷層平面中某一點(diǎn)的密度值可看作這一平面內(nèi)所有經(jīng)過(guò)該點(diǎn)的射線投影之和（的平均值）該點(diǎn)的射線投影之和（的平均值）”5.6 圖像投影重建圖像投影重建 Rado

26、n變換的matlab實(shí)現(xiàn)radon：Radon transform SyntaxR = radon(I, theta)R,xp = radon(.)Description:R = radon(I, theta) returns the Radon transform R of the intensity image I for the angle theta degrees.The Radon transform is the projection of the image intensity along a radial line oriented at a specific angle.

27、If theta is a scalar, R is a column vector containing the Radon transform for theta degrees. If theta is a vector, R is a matrix in which each column is the Radon transform for one of the angles in theta. If you omit theta, it defaults to 0:179.R,xp = radon(.) returns a vector xp containing the radi

28、al coordinates corresponding to each row of R5.6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn)% generate two imagesg1 = zeros (600,600);g1(100:500,250:350)=1;g2 = phantom (Modified Shepp-Logan,600);subplot(221);imshow(g1,);sbplot(222);imshow(g2,) % radon transformtheta = 0:0.5:179.5;R1,xp1 = radon(g1,theta);R2,xp2

29、= radon(g2,theta);5.6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn)R1 = flipud(R1); % flip up and downR2 = flipud(R2);subplot(223);imshow(R1,XData,xp1(1 end),YData,179.5 0);axis xy;axis on;xlabel(rho);ylabel(theta);subplot(224);imshow(R2,XData,xp1(1 end),YData,179.5 0);axis xy; axis on;xlabel(rho); ylabel(theta);5.

30、6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn)5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法如何利用radon變換來(lái)重建圖像f(x,y)？,( ,)( cossin,)kkkkkfx ygg xy (, ),( coscos)jjg tfx yxytdxdy 0,( , )f x yf x yd5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法第一步（first guess） =0 =2 =1 =31(0+1)5(2+3)1(0+1)5(2+3)把90角度的投影值加進(jìn)空白圖像實(shí) 例5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法第二步（second

31、 guess）021303331(0+1)5(2+3)1(0+1)5(2+3)+1(0+1)8(5+3)4(1+3)8(5+3)5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法第三步（third guess）021322443(2+1) 10(2+8)8(4+4) 12(4+8+1(0+1)8(5+3)4(1+3)8(5+35.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法第四步（fourth guess）021332133(2+1) 10(2+8)8(4+4) 12(4+8+6(3+3)12(2+10)9(1+8)15(3+12）5.6 圖像投影重建圖像投影重建 反投影重建

32、法反投影重建法6129150213所有反投影的和0/36/33/39/36-6 12-69-6 15-6063902135.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法缺點(diǎn)：星狀偽影缺點(diǎn)：星狀偽影000010000原始圖像原始圖像1/n1/n1/n1/n11/n1/n1/n1/n重構(gòu)圖像重構(gòu)圖像中心點(diǎn)中心點(diǎn)A經(jīng)經(jīng)n條投影線投影后，投條投影線投影后，投影值均為影值均為1：p1=p2=.=pn=1因此重建后因此重建后而其他點(diǎn)均為而其他點(diǎn)均為1/n：這類偽跡成為：這類偽跡成為星狀偽影星狀偽影121(.)1Anfpppn5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法缺點(diǎn)：星狀偽影

33、缺點(diǎn)：星狀偽影星狀偽影星狀偽影5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法缺點(diǎn)：圖像模糊缺點(diǎn)：圖像模糊圖像產(chǎn)生模糊5.6 圖像投影重建圖像投影重建 傅里葉切片定理（中心切片定理）傅里葉切片定理（中心切片定理）物體空間f(x,y)Radon空間 g(R)傅立葉空間F(,)RR-1F1F2F2-1中心切片定理指出：中心切片定理指出：f(x,y)f(x,y)在某一方向上的投影函數(shù)在某一方向上的投影函數(shù)g g(R)(R)的一維傅立葉變換函數(shù)的一維傅立葉變換函數(shù)G G( () )是原函數(shù)是原函數(shù)f(x,y)f(x,y)的二的二維傅立葉變換函數(shù)維傅立葉變換函數(shù)F(F(, , ) )在在( (

34、, , ) )平面上沿同一方平面上沿同一方向且過(guò)原點(diǎn)的直線上的值。向且過(guò)原點(diǎn)的直線上的值。5.6 圖像投影重建圖像投影重建 傅里葉切片定理（中心切片定理）傅里葉切片定理（中心切片定理）投影函數(shù)的數(shù)學(xué)表達(dá)式：投影函數(shù)的數(shù)學(xué)表達(dá)式： dxdyRyxyxfdlyxfRg)sincos(),(),()(f(x,y)的二維傅立葉變換：的二維傅立葉變換：)()()sincos(),(),(),(122)sincos(2RgFdReRgdxdydReRyxyxfdxdyeyxfFRjRjyxj5.6 圖像投影重建圖像投影重建 傅里葉變換法傅里葉變換法2D IFT空間域空間域頻域頻域1D FT:(0 )( ,

35、 )f x y( )G插值插值( , )F u v( , )F ( )gR5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法dudeuFyxfyuxj 2,目標(biāo)函數(shù) f(x,y) 可由傅立葉函數(shù)F(u,v) 的逆變換獲得，即5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法雅可比行列式sincosu頻域中的笛卡爾坐標(biāo)與極坐標(biāo)的關(guān)系為：sincosudddduudud 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法dudeuFyxfyuxj 2,deFddudeuFyxfyxjyuxjsincos20202sin,cos, dddudsincosu5.6 圖像投影重建

36、圖像投影重建 濾波反投影法濾波反投影法Let: F(cos, sin)=P(, ) dePddePddePddeFdyxfyxjyxjyxjyxjsincos200sincos200sincos2200sincos2200,sin,cos,P(w,) 為投影變換的一為投影變換的一維傅里葉變換維傅里葉變換5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法,tptp,PPdePddePddePdyxfyxjyxjyxjsincos20sincos200sincos200,5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法ddePyxftj20, P(, )表示對(duì)應(yīng)于角度的單位投影的

37、傅立葉變換；里層的積分是P(, )|的逆傅立葉變換，記為g(t,)，在空間域，它表示單位投影被一頻域響應(yīng)為|的函數(shù)做濾波運(yùn)算，故稱之為濾波反投影1D Fourier transforminverse 1D Fourier transformbackprojection for all anglesfilter5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法1D FT空間域空間域頻域頻域1D IFT( , )f x y1( )F gR( )gR( )gR1( )F gR濾波器濾波器5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇頻域中：濾波器|空間

38、域中與其對(duì)應(yīng)的濾波器為： dettj2將t0代入上式計(jì)算得到(0)，即曲線|以下的面積。當(dāng)時(shí)，(0)，所以上式是無(wú)法直接計(jì)算的，必須另想它法，引入限帶函數(shù)（band-limiting function）5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 2jtted濾波器是個(gè)無(wú)限頻帶的濾波函數(shù), 由于其積分是發(fā)散的, 根據(jù)佩利- - 維納準(zhǔn)則, 這一理想濾波器是不可實(shí)現(xiàn)的。實(shí)際數(shù)值計(jì)算通常采用加窗的濾波函數(shù)。運(yùn)用不同的窗函數(shù)可以得到不同的濾波器5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 10H wothers 例如

39、頻域?yàn)V波器例如頻域?yàn)V波器5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇nRam-Lak: 矩形窗nShepp-Logan正弦窗nCosine: 余弦窗nHamming: 通用Hamming窗5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇(a)圖為理想濾波器 (b)圖為修正后濾波器 亦理論上濾波器亦稱為Ramp濾波器，其高頻分量是無(wú)限延伸的，但實(shí)際實(shí)現(xiàn)時(shí)必須截?cái)嗵幚?，如圖 (b)圖中虛線所示，相當(dāng)于在帶寬之外突然衰減為零，在重建圖像的邊緣時(shí)會(huì)出現(xiàn)環(huán)狀震蕩條紋，稱之為Gibbs現(xiàn)象。為有效地消除此現(xiàn)象，我們需對(duì)Ramp

40、濾波器稍作平滑處理，如將之與作卷積，得到Shepp-Logan濾波器； 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇平滑了圖像，損失了部分高頻信息 Shepp-Logan濾波器5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇Hamming濾波器降低了高頻噪聲, 可得到Hamming濾波器和Hanning濾波器 5.6 圖像投影重建圖像投影重建反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn)I = iradon(R, theta)I = iradon(P, theta, interp, filter, frequ

41、ency_scaling, output_size)I,H = iradon(.)DescriptionI = iradon(R, theta) reconstructs the image I from projection data in the two-dimensional array R. The columns of R are parallel beam projection data. iradon assumes that the center of rotation is the center point of the projections, which is defin

42、ed as ceil(size(R,1)/2).theta describes the angles (in degrees) at which the projections were taken. It can be either a vector containing the angles or a scalar specifying D_theta, the incremental angle between projections. If theta is a vector, it must contain angles with equal spacing between them

43、. If theta is a scalar specifying D_theta, the projections were taken at angles theta = m*D_theta, where m = 0,1,2,.,size(R,2)-1. If the input is the empty matrix (), D_theta defaults to 180/size(R,2).5.6 圖像投影重建圖像投影重建反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn)I = iradon(P, theta, interp, filter, frequency_scaling, o

44、utput_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four arguments. iradon uses default values for any of these arguments that you erp specifies the type of interpolation to use in the back projection. The available options are

45、listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbor interpolation Linear: Linear interpolation (the default)spline: Spline interpolationcubic :Cubic interpolation from MATLAB 5. 5.6 圖像投影重建圖像投影重建反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn)I = iradon(P, theta, interp, filter,

46、 frequency_scaling, output_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four arguments. iradon uses default values for any of these arguments that you erp specifies the type of interpolation to use in the back projection. The a

47、vailable options are listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbor interpolation Linear: Linear interpolation (the default)spline: Spline interpolationcubic :Cubic interpolation from MATLAB 5. 5.6 圖像投影重建圖像投影重建反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn)I = iradon(P, t

48、heta, interp, filter, frequency_scaling, output_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four arguments. iradon uses default values for any of these arguments that you omit.filter specifies the filter to use for frequency domain fi

49、ltering. filter can be any of the strings that specify standard filters.Ram-LakShepp-LoganCosineHammingHannnone5.6 圖像投影重建圖像投影重建圖像重建圖像重建exampleclear;close clear;close allall; clc; clc; g = phantom (g = phantom (Modified Shepp-LoganModified Shepp-Logan,600);,600);subplot(241);imshow(g,);title (subplot

50、(241);imshow(g,);title (original imageoriginal image) )theta = 0:0.5:179.5;theta = 0:0.5:179.5;R,xp = radon(g,theta);R,xp = radon(g,theta);subplot(242);imshow(R,);title (subplot(242);imshow(R,);title (Projection imageProjection image) )f1 = iradon(R,theta,f1 = iradon(R,theta,nonenone); );subplot(243

51、);imshow(f1,);title (subplot(243);imshow(f1,);title (Reconstructed Image with back Reconstructed Image with back projectionprojection) )f2 = iradon(R,theta,f2 = iradon(R,theta,Ram-LakRam-Lak); );subplot(244);imshow(f2,);title (subplot(244);imshow(f2,);title (Reconstructed Image with Ram-Reconstructe

52、d Image with Ram-Lak filterLak filter) )5.6 圖像投影重建圖像投影重建圖像重建圖像重建examplef3 = iradon(R,theta,f3 = iradon(R,theta,Shepp-LoganShepp-Logan); );subplot(245);imshow(f3,);title (subplot(245);imshow(f3,);title (Reconstructed Image Reconstructed Image with Shepp-Logan filterwith Shepp-Logan filter) )f4 = irad

53、on(R,theta,f4 = iradon(R,theta,CosineCosine); );subplot(246);imshow(f4,);title (subplot(246);imshow(f4,);title (Reconstructed Image Reconstructed Image with Cosine filterwith Cosine filter) )f5 = iradon(R,theta,f5 = iradon(R,theta,HammingHamming); );subplot(247);imshow(f5,);title (subplot(247);imsho

54、w(f5,);title (Reconstructed Image Reconstructed Image with Hamming filterwith Hamming filter) )f6 = iradon(R,theta,f6 = iradon(R,theta,HannHann); );subplot(248);imshow(f6,);title (subplot(248);imshow(f6,);title (Reconstructed Image Reconstructed Image with Hann filterwith Hann filter) )5.6 圖像投影重建圖像投

55、影重建5.6 圖像投影重建圖像投影重建扇形掃描重建扇形掃描重建成像幾何5.6 圖像投影重建圖像投影重建扇形掃描重建扇形掃描重建扇束情況下的重建算扇束情況下的重建算法較為復(fù)雜，但實(shí)質(zhì)法較為復(fù)雜，但實(shí)質(zhì)沒(méi)有改變?？刹捎闷?jīng)]有改變?？刹捎闷叫惺闆r下的算法實(shí)行束情況下的算法實(shí)現(xiàn)，只需加以適當(dāng)?shù)噩F(xiàn)，只需加以適當(dāng)?shù)匦拚纯尚拚纯蒼 重排算法重排算法: : 把一個(gè)視圖中采得的扇形數(shù)據(jù)重新組合成平行的把一個(gè)視圖中采得的扇形數(shù)據(jù)重新組合成平行的 射線投影數(shù)據(jù)，然后采用平行束重建算法重建射線投影數(shù)據(jù)，然后采用平行束重建算法重建n 直接重建算法直接重建算法: : 不必?cái)?shù)據(jù)重排，只需適當(dāng)加權(quán)即可運(yùn)用與平不必?cái)?shù)據(jù)重

56、排，只需適當(dāng)加權(quán)即可運(yùn)用與平 行束類似的算法重建行束類似的算法重建5.6 圖像投影重建圖像投影重建扇形束重建- EA等角扇形束重建n當(dāng)同樣大小的探測(cè)器單元沿著中心為X射線焦點(diǎn)的弧排列時(shí)，就形成等角采樣n扇形束的每一條射線可由和確定，其中是射線與中心射線（假想的通過(guò)X射線源和等中心的直線）的夾角，稱為探測(cè)器角；是中心射線與y軸的夾角，稱為投影角 5.6 圖像投影重建圖像投影重建扇形束重建- EA等角扇形束重建n 投影乘以探測(cè)器角的余弦，濾波后的樣本隨著到光源的距離的增長(zhǎng)而增長(zhǎng)n 重建公式可由用(t, ) 坐標(biāo)確定(, ) 坐標(biāo)上的每個(gè)樣本來(lái)得到。扇形投影中的投影樣本q(, )就轉(zhuǎn)化為平行投影中的

57、投影樣本p(t, ) 5.6 圖像投影重建圖像投影重建扇形投影和重建的扇形投影和重建的matlab實(shí)現(xiàn)實(shí)現(xiàn)fanbeam：Fan-beam transform Syntax:F = fanbeam(I,D)F = fanbeam(., param1, val1, param1, val2,.)F, fan_sensor_positions, fan_rotation_angles = fanbeam(.)FanRotationIncrement - Positive real scalar specifying the increment of the rotation angle of th

58、e fan-beam projections. Measured in degrees. Default value is 1.FanSensorGeometry - Text string specifying how sensors are positioned. Valid values are arc or line. In the arc geometry, sensors are spaced equally along a circular arc, This is the default value.FanSensorSpacing - Positive real scalar

59、 specifying the spacing of the fan-beam sensors. Interpretation of the value depends on the setting of FanSensorGeometry. If FanSensorGeometry is set to arc (the default), the value defines the angular spacing in degrees. Default value is 1. If FanSensorGeometry is line, the value specifies the line

60、ar spacing. Default value is 1.5.6 圖像投影重建圖像投影重建扇形投影和重建的扇形投影和重建的matlab實(shí)現(xiàn)實(shí)現(xiàn)fanbeam：Fan-beam transform Syntax:F = fanbeam(I,D)F = fanbeam(., param1, val1, param1, val2,.)F, fan_sensor_positions, fan_rotation_angles = fanbeam(.)FanRotationIncrement - Positive real scalar specifying the increment of the