第五章 圖像復(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 圖像投影重建圖像投影重建 周期噪聲的模型是二維正弦波，通過帶阻、帶通和陷波濾波器 可以被有效去除。 理想帶阻濾波器的表達(dá)式理想帶阻濾波器的表達(dá)式: : 0 00 0 1,( , ) 2 ( , )0,( , ) 22 1,( , ) 2 W D

2、 u vD WW H u vDD u vD W D u vD 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 帶阻濾波器帶阻濾波器 n n階的巴特沃思帶阻濾波器階的巴特沃思帶阻濾波器 2 22 0 1 ( , ) ( , ) 1 ( , ) n H u v D u v W D u vD 高斯帶阻濾波器高斯帶阻濾波器 2 22 0 ( , )1 2( , ) ( , )1 Du vD D u v W H u ve 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 帶阻濾波器帶阻濾波器 (a)理想帶阻濾波器理想帶阻濾波器 (b)巴特沃思帶阻濾波器巴特沃思帶阻濾波器 (c)高斯帶阻濾波器高斯帶阻濾

3、波器 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 帶阻濾波器帶阻濾波器 (a) 被正弦噪聲污染的圖像被正弦噪聲污染的圖像 (b) 圖圖(a)的頻譜的頻譜 (c) 巴特沃思帶阻濾波器巴特沃思帶阻濾波器 (d) 濾波效果圖濾波效果圖 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 ( , )1( , ) bpbr Hu vHu v 帶通濾波器帶通濾波器 帶通濾波器執(zhí)行與帶阻濾波器相反的操作。帶通濾波器執(zhí)行與帶阻濾波器相反的操作。 ( , ) ( , ): bp br Hu v Hu v 帶通濾波器的傳遞函數(shù)可根據(jù)相應(yīng)的 帶阻濾波器的傳遞函數(shù)得到 不直接使用，損失大量不直接使用，損失大量 圖

4、像細(xì)節(jié)。圖像細(xì)節(jié)。 可利用帶通濾波器提取噪可利用帶通濾波器提取噪 聲模式。聲模式。 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 陷波濾波器陷波濾波器 阻止阻止(或通過或通過)事先定義的中心頻率鄰域內(nèi)的頻率。事先定義的中心頻率鄰域內(nèi)的頻率。 (a) 理想陷波濾波器理想陷波濾波器 (b) 巴特沃思陷波濾波器巴特沃思陷波濾波器 (c) 高斯陷波濾波器高斯陷波濾波器 由于傅立葉變換由于傅立葉變換 是對(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 ,(

5、,)(,)Du vuv半徑為中心在且在對(duì)稱的理想陷波濾波器的傳遞函數(shù) 1020 0( , )( , ) ( , ) 1 D u vDD u vD H u v 或 其他 22 1/2 100 22 1/2 200 ( , )(/2)(/2) ( , )(/2)(/2) D u vuMuvNv D u vuMuvNv 其中 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 陷波濾波器陷波濾波器 2 0 12 : 1 ( , ) 1 ( , )/, n n H u v D D u vDu v 階數(shù)為 的巴特沃思陷波帶阻濾波器的傳遞函數(shù)為 12 2 0 ( , )/,1 2 : ( , )1 D u

6、vDu v D H u ve 高斯陷波帶阻濾波器的傳遞函數(shù)為 還可以得到另一種陷波濾波器還可以得到另一種陷波濾波器,它能通過它能通過 (而不是阻止而不是阻止)包含在陷波區(qū)的頻率包含在陷波區(qū)的頻率. 陷波區(qū)域的形狀可以是任意的陷波區(qū)域的形狀可以是任意的(如矩形如矩形)。 5.4 頻域?yàn)V波降低周期噪聲頻域?yàn)V波降低周期噪聲 圖像退化模型：圖像退化模型： 5.5 退化函數(shù)建模退化函數(shù)建模 ),(),(*),(),(yxyxfyxhyxg ),(),(),(),(vuNvuFvuHvuG 退化系統(tǒng)一般情況下是：線性，位置不變的退化系統(tǒng)退化系統(tǒng)一般情況下是：線性，位置不變的退化系統(tǒng) （1 1）線性：）線性

7、： yxfbHyxfaHyxbfyxafH, 2121 （2 2）位置不變性：對(duì)任意）位置不變性：對(duì)任意 , yxf 有有 yxgyxfH, 對(duì)于線性位置不變退化，圖像復(fù)原其實(shí)就是一個(gè)圖像反卷積對(duì)于線性位置不變退化，圖像復(fù)原其實(shí)就是一個(gè)圖像反卷積 的過程的過程 圖像觀察估計(jì)法圖像觀察估計(jì)法 給定一幅退化圖像，但沒有退化函數(shù) H 的知識(shí)，那么估計(jì) 該函數(shù)的方法之一就是收集圖像自身的信息： 尋找簡(jiǎn)單結(jié)構(gòu)的子圖像 尋找受噪聲影響小的子圖像 5.5 退化函數(shù)建模退化函數(shù)建模 估計(jì)退化系統(tǒng)模型的三種方法估計(jì)退化系統(tǒng)模型的三種方法 構(gòu)造一個(gè)估計(jì)圖像，它與觀察的子圖像有相同大小和特性 表示觀察子圖像， 表示

8、構(gòu)造的子圖像 和 為對(duì)應(yīng)的傅立葉變換。 ( , )( , )?/( , ) S ss H u vGu vF u v ),(yxgs),( yxfs ),(vuGs),( vuFs 假設(shè)空間不變的，由 推導(dǎo)出完全函數(shù) ),(vuH s ),(vuH 5.5 退化函數(shù)建模退化函數(shù)建模 圖像試驗(yàn)估計(jì)法圖像試驗(yàn)估計(jì)法 使用與被退化圖像設(shè)備相似的裝置，并得到一個(gè)脈沖的沖激 響應(yīng)，可以進(jìn)行較準(zhǔn)確的退化估計(jì)： ( , ) ( , ) G u v H u v A 一個(gè)脈沖點(diǎn)一個(gè)脈沖點(diǎn)成像系統(tǒng)成像系統(tǒng)H 此處A是一個(gè)沖激的傅立葉變換，表示沖擊強(qiáng)度，為一常數(shù)。 右圖為一個(gè)放大的 亮脈沖以及退化的 沖激。 ( ,

9、)g x y 退化圖像退化圖像 模型估計(jì)法模型估計(jì)法 建立退化模型，考慮引起退化的環(huán)境因素。建立退化模型，考慮引起退化的環(huán)境因素。 22 5/6 () ( , ) k uv H u ve 例如：例如：Hufnagel Hufnagel 等等 Stanley Stanley 的退化模型就是基于大氣湍的退化模型就是基于大氣湍 流的物理特性而提出來的，其中流的物理特性而提出來的，其中k k為常數(shù)，與湍流特性相關(guān)。為常數(shù)，與湍流特性相關(guān)。 ( (除了指數(shù)除了指數(shù)5/65/6，該公式與高斯低通濾波形式相同，該公式與高斯低通濾波形式相同.).) 5.5 退化函數(shù)建模退化函數(shù)建模 模型估計(jì)法模型估計(jì)法 5.

10、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 uv H u ve 如果已知系統(tǒng)的傳遞函數(shù) ，則根據(jù) vuFvuHvuG, vuH, vuHvuGvuF, 可得復(fù)原圖像的譜，經(jīng)傅氏逆變換即可得到復(fù)原圖像 在忽略噪聲的影響，退化模型的傅氏變換為 實(shí)際應(yīng)用時(shí)存在病態(tài)的問題，即在 H(u,v) 等于零或非常小

11、 的數(shù)值點(diǎn)上， 將變成無窮大或非常大的數(shù)。 ),( vuF - 這就是逆濾波復(fù)原法 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波 vuNvuFvuHvuG, vuH vuN vuF vuH vuN vuH vuG vuF , , ),( , , , , , 系統(tǒng)中存在噪聲時(shí)退化模型的傅立葉變換為： 寫成逆濾波復(fù)原的方式： 1）即使知道退化函數(shù)，也不能準(zhǔn)確復(fù)原圖像，因?yàn)樵肼暫瘮?shù) N(u,v) 是 一個(gè)隨機(jī)函數(shù)，其傅里葉變換未知。 2）如果退化是零或非常小的值，噪聲即使數(shù)值很小，但 N(u,v)/H(u,v) 之比 (上式第二項(xiàng)) 可能非常大，很容易錯(cuò)估 的值。),( vuF 12 () (

12、 , )( , )( , )( , ) jux vy f 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) 行，即將退化圖像的傅立葉譜限制在沒出現(xiàn)零點(diǎn)而 且數(shù)值又不是太小的有限范圍內(nèi)，即通過將頻率限 制為接近原點(diǎn)分析，減少了遇到零值的幾率。 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-逆濾波逆濾波 劇烈湍流劇烈湍流( (k k=0.0025)=0.0025) 大氣湍流模型模擬退化模糊一幅圖

13、像大氣湍流模型模擬退化模糊一幅圖像 可忽略的湍流可忽略的湍流 6522 22 / )/()/( ),( NvMuk evuH 對(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é)果，一般直 接逆濾波效果較差。接逆濾波效果

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

15、考慮退化函數(shù)和噪聲統(tǒng)計(jì)特征。 E期望值。期望值。 ),(yxf 因此維納濾波復(fù)原又稱為最小均方誤差復(fù)原。因此維納濾波復(fù)原又稱為最小均方誤差復(fù)原。 min),(),( 2 yxfyxfE 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原 ),( ),(/ ),(),( ),( ),( ),( vuG vuSvuSvuH vuH vuH vuF f 2 2 1 f 誤差函數(shù)的最小值在頻域里可以通過近似圖像 的 傅里葉變換來計(jì)算： 葉變換為復(fù)原近似圖像的傅里 換為退化圖像的傅里葉變 換為退化函數(shù)的傅里葉變其中： ),( ),( ),( vuF vuG vuH 為未退化圖像的功率譜

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

17、定的常數(shù)。是根據(jù)信噪比的某種先驗(yàn)知識(shí)確定的常數(shù)。 ),(/ ),(),( ),( ),( ),( vuSvuSvuH vuH vuH vuH f w 2 2 1 維納濾波復(fù)原：維納濾波復(fù)原： 維納濾波需要假定下述條件成立維納濾波需要假定下述條件成立( (或近似成立或近似成立) )： (1)(1) 系統(tǒng)為線性、空間不變；系統(tǒng)為線性、空間不變； (2)(2) 退化圖像、原始圖像和噪聲都是均勻隨機(jī)場(chǎng)，噪聲的均退化圖像、原始圖像和噪聲都是均勻隨機(jī)場(chǎng)，噪聲的均 值為零，且與圖像不相關(guān)。值為零，且與圖像不相關(guān)。 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原 維納濾波復(fù)原與逆濾波復(fù)原的

18、比較 全頻逆濾波半徑受限逆濾波維納濾波復(fù)原 (交互選擇K) 維納濾波的缺點(diǎn)： 未退化圖像和噪聲的功率譜必須是已知的； 功率比(信噪比)常數(shù)K 的估計(jì)一般還是沒有合適的解。 5.6 圖像復(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 Syntax J = deconvwnr(I,PSF) J = deconvwnr(I,PSF,NSR) J = deconvwnr(I,PSF,NCORR,

19、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ù)原的方法-Wiener濾波復(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(

20、f,); title(title(orginal clean imageorginal clean image); ); % generate the degrade function% generate the degrade function PSF = fspecial(PSF = fspecial(motionmotion,7,45);,7,45); subplot(232); imshow(PSF,);subplot(232); imshow(PSF,); title(title(Point spread functionPoint spread function); ); % us

21、ing the PSF to degrade image% using the PSF to degrade image gb = imfilter(f,PSF,gb = imfilter(f,PSF,circularcircular); ); 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ù)原案例

22、：維納濾波和逆濾波復(fù)原案例： % add noise to the degraded image% add noise to the degraded image noise = imnoise(zeros(size(f),noise = 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 ima

23、ge with noise) ) % inverse filtering% inverse filtering fr1 = deconvwnr(g,PSF);fr1 = deconvwnr(g,PSF); subplot(235);imshow(fr1,)subplot(235);imshow(fr1,) title(title(inverse filtering resultinverse filtering result) ) 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原 維納濾波和逆濾波復(fù)原案例：維納濾波和逆濾波復(fù)原案例： % wiener filtering% w

24、iener filtering Sn = abs(fft2(noise).2; % noise power spectrum nA = sum(Sn(:)/prod(size(noise); % noise average power Sf = abs(fft2(f).2; % image power spectrum fA = sum(Sf(:)/prod(size(f); % image average power R = nA/fA; % signal to noise ratio fr2 = deconvwnr(g,PSF,R);fr2 = deconvwnr(g,PSF,R); su

25、bplot(236);imshow(fr2,)subplot(236);imshow(fr2,) title(title(wiener filtering resultwiener filtering result) ) 5.6 圖像復(fù)原的方法圖像復(fù)原的方法-Wiener濾波復(fù)原濾波復(fù)原 維納濾波和逆濾波復(fù)原案例：維納濾波和逆濾波復(fù)原案例： 5.6 圖像投影重建圖像投影重建 概念：投影重建一般指利用物體的多個(gè)（軸向概念：投影重建一般指利用物體的多個(gè)（軸向 ）投影圖像重建目標(biāo)圖像的過程。它是一類特）投影圖像重建目標(biāo)圖像的過程。它是一類特 殊的圖像處理方法，輸入的是一系列的投影圖殊的圖像處理方法，

26、輸入的是一系列的投影圖 ，輸出是重建圖。，輸出是重建圖。 通過投影重建可以直接的看到原來被投影的物通過投影重建可以直接的看到原來被投影的物 體的某種特性的空間分布，比直觀觀測(cè)投影圖體的某種特性的空間分布，比直觀觀測(cè)投影圖 要直觀的多。要直觀的多。 5.6 圖像投影重建圖像投影重建 Radon變換 對(duì)f(x,y)的Radon變換g(t, )定義為沿由t和 定義的直線l的線積分。 ,( , ) ,( cossin) kkj g tf x ydl f x yxyt dxdy 5.6 圖像投影重建圖像投影重建 Radon變換 (, ),( coscos) jj g tfx yxytdxdy Radon

27、 變換揭示了函數(shù)和投影之間的關(guān)系，若函數(shù)為f (x, y)， 則不同角度下的投影可寫為 原理：原理：“斷層平面中某一點(diǎn)的密度值可看作這一平面內(nèi)所有經(jīng)過斷層平面中某一點(diǎn)的密度值可看作這一平面內(nèi)所有經(jīng)過 該點(diǎn)的射線投影之和（的平均值）該點(diǎn)的射線投影之和（的平均值）” 5.6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn) radon：Radon transform Syntax R = radon(I, theta) R,xp = radon(.) Description: R = radon(I, theta) returns the Radon transform R of the

28、 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. 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 eac

29、h 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 radial coordinates corresponding to each row of R 5.6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn) % generate two images g1 = zeros (600,600); g1(100:500,250:3

30、50)=1; g2 = phantom (Modified Shepp-Logan,600); subplot(221);imshow(g1,); sbplot(222);imshow(g2,) % radon transform theta = 0:0.5:179.5; R1,xp1 = radon(g1,theta); R2,xp2 = radon(g2,theta); 5.6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn) R1 = flipud(R1); % flip up and down R2 = flipud(R2); subplot(223); imshow(R1,

31、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.6 圖像投影重建圖像投影重建 Radon變換的matlab實(shí)現(xiàn) 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 如何利用radon變換來重建圖像f(x,y)？ ,( ,)( cossin ,) k kkkk fx ygg

32、 xy (, ),( coscos) jj g tfx yxytdxdy 0 ,( , )f x yf x yd 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 第一步（first guess） =0 =2 =1 =3 1(0+1)5(2+3) 1(0+1)5(2+3) 把90角度的投影值加進(jìn)空白圖像 實(shí) 例 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 第二步（second guess） 02 13 03 33 1(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 圖像投影重建圖像投影重建 反投影重建法反

33、投影重建法 第三步（third guess） 02 13 22 44 3(2+1) 10(2+8) 8(4+4) 12(4+8 + 1(0+1)8(5+3) 4(1+3)8(5+3 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 第四步（fourth guess） 02 13 32 13 3(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 圖像投影重建圖像投影重建 反投影重建法反投影重建法 612 915 02 13 所有反投影的和 0/36/3 3/39/3 6-6 12-6 9-6 15-6 06 3

34、9 02 13 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 缺點(diǎn)：星狀偽影缺點(diǎn)：星狀偽影 000 010 000 原始圖像原始圖像 1/n1/n1/n 1/n11/n 1/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：這類偽跡成為：這類偽跡成為 星狀偽影星狀偽影 12 1 (.)1 An fppp n 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 缺點(diǎn)：星狀偽影缺點(diǎn)：星狀偽影 星狀偽影星狀偽影 5.6 圖像投影重建圖像

35、投影重建 反投影重建法反投影重建法 缺點(diǎn)：圖像模糊缺點(diǎn)：圖像模糊 圖像產(chǎn)生模糊 5.6 圖像投影重建圖像投影重建 傅里葉切片定理（中心切片定理）傅里葉切片定理（中心切片定理） 物體空間f(x,y) Radon空間 g(R) 傅立葉空間F(,) R R-1 F1 F2 F2-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(, , ) )在在( (, , ) )平面

36、上沿同一方平面上沿同一方 向且過原點(diǎn)的直線上的值。向且過原點(diǎn)的直線上的值。 5.6 圖像投影重建圖像投影重建 傅里葉切片定理（中心切片定理）傅里葉切片定理（中心切片定理） 投影函數(shù)的數(shù)學(xué)表達(dá)式：投影函數(shù)的數(shù)學(xué)表達(dá)式： dxdyRyxyxf dlyxfRg )sincos(),( ),()( f(x,y)的二維傅立葉變換：的二維傅立葉變換： )( )( )sincos(),( ),(),( 1 2 2 )sincos(2 RgF dReRg dxdydReRyxyxf dxdyeyxfF Rj Rj yxj 5.6 圖像投影重建圖像投影重建 傅里葉變換法傅里葉變換法 2D IFT 空間域空間域

37、頻域頻域 1D FT :(0 ) ( , )f x y ( )G 插值插值 ( , )F u v ( , )F ( )gR 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 dudeuFyxf yuxj 2 , 目標(biāo)函數(shù) f(x,y) 可由傅立葉函數(shù)F(u,v) 的逆 變換獲得，即 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 雅可比行列式 sin cos u 頻域中的笛卡爾坐標(biāo)與 極坐標(biāo)的關(guān)系為： sin cos u dddd uu dud 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 dudeuFyxf yuxj 2 , deFd dudeuFyxf

38、yxj yuxj sincos2 0 2 0 2 sin,cos , dddud sin cos u 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 Let: F(cos, sin)=P(, ) dePd dePd dePd deFdyxf yxj yxj yxj yxj sincos2 00 sincos2 00 sincos2 2 00 sincos2 2 00 , , , sin,cos, P(w,) 為投影變換的一為投影變換的一 維傅里葉變換維傅里葉變換 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 ,tptp ,PP dePd dePd dePdyxf yx

39、j yxj yxj sincos2 0 sincos2 00 sincos2 00 , , , 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 ddePyxf tj2 0 , P(, )表示對(duì)應(yīng)于角度的單位投影的傅立葉變換；里層的積分是 P(, )|的逆傅立葉變換，記為g(t,)，在空間域，它表示單位 投影被一頻域響應(yīng)為|的函數(shù)做濾波運(yùn)算，故稱之為濾波反投 影 1D Fourier transform inverse 1D Fourier transform backprojection for all angles filter 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波

40、反投影法 1D FT 空間域空間域 頻域頻域 1D IFT ( , )f x y 1 ( )F gR ( )gR ( )gR 1 ( )F gR 濾波器濾波器 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 頻域中：濾波器| 空間域中與其對(duì)應(yīng)的濾波器為： det tj 2 將t0代入上式計(jì)算得到(0)，即曲線|以下的面積。 當(dāng)時(shí)，(0)，所以上式是無法直接計(jì)算的，必須 另想它法，引入限帶函數(shù)（band-limiting function） 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 2jt ted 濾波器是個(gè)無

41、限頻帶的濾波函數(shù), 由于其積分是發(fā)散的, 根 據(jù)佩利- - 維納準(zhǔn)則, 這一理想濾波器是不可實(shí)現(xiàn)的。實(shí)際 數(shù)值計(jì)算通常采用加窗的濾波函數(shù)。運(yùn)用不同的窗函數(shù)可以 得到不同的濾波器 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 1 0 H w others 例如頻域?yàn)V波器例如頻域?yàn)V波器 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 nRam-Lak: 矩形窗 nShepp-Logan正弦窗 nCosine: 余弦窗 nHamming: 通用Hamming窗 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選

42、擇濾波反投影法中濾波器的選擇 (a)圖為理想濾波器 (b)圖為修正后濾波器 亦 理論上濾波器亦稱為Ramp濾波器，其高頻分量是無限延伸的，但 實(shí)際實(shí)現(xiàn)時(shí)必須截?cái)嗵幚?，如圖 (b)圖中虛線所示，相當(dāng)于在帶 寬之外突然衰減為零，在重建圖像的邊緣時(shí)會(huì)出現(xiàn)環(huán)狀震蕩條紋， 稱之為Gibbs現(xiàn)象。為有效地消除此現(xiàn)象，我們需對(duì)Ramp濾波器稍 作平滑處理，如將之與作卷積，得到Shepp-Logan濾波器； 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 平滑了圖像，損失了部分高頻信息 Shepp-Logan濾波器 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器

43、的選擇濾波反投影法中濾波器的選擇 Hamming濾波器 降低了高頻噪聲, 可得到Hamming濾波 器和Hanning濾波器 5.6 圖像投影重建圖像投影重建 反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn) I = iradon(R, theta) I = iradon(P, theta, interp, filter, frequency_scaling, output_size) I,H = iradon(.) Description I = iradon(R, theta) reconstructs the image I from projection data in the t

44、wo- 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 defined as ceil(size(R,1)/2). theta describes the angles (in degrees) at which the projections were taken. It can be either a vecto

45、r 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. 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

46、)- 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, output_size) specifies parameters to use in the inverse Radon transform. You can specify any combination of the last four

47、arguments. iradon uses default values for any of these arguments that you omit. interp specifies the type of interpolation to use in the back projection. The available options are listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbor interpolation Linear: Li

48、near interpolation (the default) spline: Spline interpolation cubic :Cubic interpolation from MATLAB 5. 5.6 圖像投影重建圖像投影重建 反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn) I = iradon(P, theta, interp, filter, frequency_scaling, output_size) specifies parameters to use in the inverse Radon transform. You can specify any com

49、bination of the last four arguments. iradon uses default values for any of these arguments that you omit. interp specifies the type of interpolation to use in the back projection. The available options are listed in order of increasing accuracy and computational complexity. nearest : Nearest-neighbo

50、r interpolation Linear: Linear interpolation (the default) spline: Spline interpolation cubic :Cubic interpolation from MATLAB 5. 5.6 圖像投影重建圖像投影重建 反投影重建算法的反投影重建算法的matlab實(shí)現(xiàn)實(shí)現(xiàn) I = iradon(P, theta, interp, filter, frequency_scaling, output_size) specifies parameters to use in the inverse Radon transfor

51、m. 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 filtering. filter can be any of the strings that specify standard filters. Ram-Lak Shepp-Logan Cosine Hamming Ha

52、nn none 5.6 圖像投影重建圖像投影重建 圖像重建圖像重建example clear;close clear;close allall; clc; clc; g = phantom (g = phantom (Modified Shepp-LoganModified Shepp-Logan,600);,600); subplot(241);imshow(g,);title (subplot(241);imshow(g,);title (original imageoriginal image) ) theta = 0:0.5:179.5;theta = 0:0.5:179.5; R,x

53、p = 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);imshow(f1,);title (subplot(243);imshow(f1,);title (Reconstructed Image with back Reconstructe

54、d 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-Reconstructed Image with Ram- Lak filterLak filter) ) 5.6 圖像投影重建圖像投影重建 圖像重建圖像重建example f3 = iradon(R,thet

55、a,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 = iradon(R,theta,f4 = iradon(R,theta,CosineCosine); ); subplot(246);imshow(f4,);title (subplo

56、t(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);imshow(f5,);title (Reconstructed Image Reconstructed Image with Hamming filterwith Hammin

57、g 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 圖像投影重建圖像投影重建 5.6 圖像投影重建圖像投影重建 扇形掃描重建扇形掃描重建 成像幾何 5.6 圖像投影重建圖像投影重建 扇形掃描重建扇形掃描重建 扇束情況下的重建算扇束情

58、況下的重建算 法較為復(fù)雜，但實(shí)質(zhì)法較為復(fù)雜，但實(shí)質(zhì) 沒有改變?？刹捎闷?jīng)]有改變?？刹捎闷?行束情況下的算法實(shí)行束情況下的算法實(shí) 現(xiàn)，只需加以適當(dāng)?shù)噩F(xiàn)，只需加以適當(dāng)?shù)?修正即可修正即可 n 重排算法重排算法: : 把一個(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ù)重排，只需適當(dāng)加權(quán)即可運(yùn)用與平 行束類似的算法重建行束類似的算法重建 5.6 圖像投影重建圖像投影重建 扇形束重建- EA等角扇形束重建

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

60、b實(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 the fan-beam projections. Measured in degrees