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

1、數(shù)字圖像處理數(shù)字圖像處理 主講人：杜宏博 第五章第五章 圖像復原與重建圖像復原與重建 5.1 圖像退化圖像退化/復原處理的模型復原處理的模型 5.2 噪聲模型噪聲模型 5.3 空間濾波去噪空間濾波去噪 5.4 頻域濾波去噪頻域濾波去噪 5.5 退化函數(shù)建模退化函數(shù)建模 5.6 圖像復原的方法圖像復原的方法 直接逆濾波直接逆濾波 維納濾波維納濾波 5.7 圖像投影重建圖像投影重建 周期噪聲的模型是二維正弦波，通過帶阻、帶通和陷波濾波器 可以被有效去除。 理想帶阻濾波器的表達式理想帶阻濾波器的表達式: : 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 頻域濾波降低周期噪聲頻域濾波降低周期噪聲 帶阻濾波器帶阻濾波器 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 頻域濾波降低周期噪聲頻域濾波降低周期噪聲 帶阻濾波器帶阻濾波器 (a)理想帶阻濾波器理想帶阻濾波器 (b)巴特沃思帶阻濾波器巴特沃思帶阻濾波器 (c)高斯帶阻濾波器高斯帶阻濾

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

4、像細節(jié)。圖像細節(jié)。 可利用帶通濾波器提取噪可利用帶通濾波器提取噪 聲模式。聲模式。 5.4 頻域濾波降低周期噪聲頻域濾波降低周期噪聲 陷波濾波器陷波濾波器 阻止阻止(或通過或通過)事先定義的中心頻率鄰域內(nèi)的頻率。事先定義的中心頻率鄰域內(nèi)的頻率。 (a) 理想陷波濾波器理想陷波濾波器 (b) 巴特沃思陷波濾波器巴特沃思陷波濾波器 (c) 高斯陷波濾波器高斯陷波濾波器 由于傅立葉變換由于傅立葉變換 是對稱的是對稱的,因此因此 陷波濾波器必須陷波濾波器必須 以關于原點對稱以關于原點對稱 的形式出現(xiàn)。的形式出現(xiàn)。 5.4 頻域濾波降低周期噪聲頻域濾波降低周期噪聲 陷波濾波器陷波濾波器 00000 ,(

5、,)(,)Du vuv半徑為中心在且在對稱的理想陷波濾波器的傳遞函數(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 頻域濾波降低周期噪聲頻域濾波降低周期噪聲 陷波濾波器陷波濾波器 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 頻域濾波降低周期噪聲頻域濾波降低周期噪聲 圖像退化模型：圖像退化模型： 5.5 退化函數(shù)建模退化函數(shù)建模 ),(),(*),(),(yxyxfyxhyxg ),(),(),(),(vuNvuFvuHvuG 退化系統(tǒng)一般情況下是：線性，位置不變的退化系統(tǒng)退化系統(tǒng)一般情況下是：線性，位置不變的退化系統(tǒng) （1 1）線性：）線性

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

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

9、)g x y 退化圖像退化圖像 模型估計法模型估計法 建立退化模型，考慮引起退化的環(huán)境因素。建立退化模型，考慮引起退化的環(huán)境因素。 22 5/6 () ( , ) k uv H u ve 例如：例如：Hufnagel Hufnagel 等等 Stanley Stanley 的退化模型就是基于大氣湍的退化模型就是基于大氣湍 流的物理特性而提出來的，其中流的物理特性而提出來的，其中k k為常數(shù)，與湍流特性相關。為常數(shù)，與湍流特性相關。 ( (除了指數(shù)除了指數(shù)5/65/6，該公式與高斯低通濾波形式相同，該公式與高斯低通濾波形式相同.).) 5.5 退化函數(shù)建模退化函數(shù)建模 模型估計法模型估計法 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, 可得復原圖像的譜，經(jīng)傅氏逆變換即可得到復原圖像 在忽略噪聲的影響，退化模型的傅氏變換為 實際應用時存在病態(tài)的問題，即在 H(u,v) 等于零或非常小

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

12、 , )( , )( , )( , ) jux vy f x yf x yN u v Hu vedudv 5.6 圖像復原的方法圖像復原的方法-逆濾波逆濾波 解決退化是零或非常小的值的途徑： 限制濾波的頻率，使其接近原點值。 在頻率平面離原點較遠的地方，H(u,v)數(shù)值較小 或為零，因此圖像復原在原點周圍的有限區(qū)域內(nèi)進 行，即將退化圖像的傅立葉譜限制在沒出現(xiàn)零點而 且數(shù)值又不是太小的有限范圍內(nèi)，即通過將頻率限 制為接近原點分析，減少了遇到零值的幾率。 5.6 圖像復原的方法圖像復原的方法-逆濾波逆濾波 劇烈湍流劇烈湍流( (k k=0.0025)=0.0025) 大氣湍流模型模擬退化模糊一幅圖

13、像大氣湍流模型模擬退化模糊一幅圖像 可忽略的湍流可忽略的湍流 6522 22 / )/()/( ),( NvMuk evuH 對退化函數(shù)對退化函數(shù)H H( (u u, ,v v) )進行精確取反并進行逆濾波，結(jié)果如下圖。進行精確取反并進行逆濾波，結(jié)果如下圖。 5.6 圖像復原的方法圖像復原的方法-逆濾波逆濾波 全頻直接全頻直接 逆濾波復原逆濾波復原 半徑為半徑為 4040時截止時截止H H 半徑為半徑為 7070時截止時截止H H 半徑為半徑為 8585時截止時截止H H 結(jié)果表明：噪聲明顯影響結(jié)果表明：噪聲明顯影響 了圖像復原結(jié)果，一般直了圖像復原結(jié)果，一般直 接逆濾波效果較差。接逆濾波效果

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

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

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

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

18、比較 全頻逆濾波半徑受限逆濾波維納濾波復原 (交互選擇K) 維納濾波的缺點： 未退化圖像和噪聲的功率譜必須是已知的； 功率比(信噪比)常數(shù)K 的估計一般還是沒有合適的解。 5.6 圖像復原的方法圖像復原的方法-Wiener濾波復原濾波復原 5.6 圖像復原的方法圖像復原的方法-Wiener濾波復原濾波復原 維納濾波器的維納濾波器的matlab實現(xiàn)實現(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ù)（點擴散函數(shù)） NSR信噪比 NCORR：噪聲的自相關函數(shù) ICORR：退化圖像的自相關函數(shù) J：復原圖像 5.6 圖像復原的方法圖像復原的方法-Wiener濾波復原濾波復原 維納濾波和逆濾波復原案例：維納濾波和逆濾波復原案例： 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 圖像復原的方法圖像復原的方法-Wiener濾波復原濾波復原 維納濾波和逆濾波復原案例

22、：維納濾波和逆濾波復原案例： % 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 圖像復原的方法圖像復原的方法-Wiener濾波復原濾波復原 維納濾波和逆濾波復原案例：維納濾波和逆濾波復原案例： % 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 圖像復原的方法圖像復原的方法-Wiener濾波復原濾波復原 維納濾波和逆濾波復原案例：維納濾波和逆濾波復原案例： 5.6 圖像投影重建圖像投影重建 概念：投影重建一般指利用物體的多個（軸向概念：投影重建一般指利用物體的多個（軸向 ）投影圖像重建目標圖像的過程。它是一類特）投影圖像重建目標圖像的過程。它是一類特 殊的圖像處理方法，輸入的是一系列的投影圖殊的圖像處理方法，

26、輸入的是一系列的投影圖 ，輸出是重建圖。，輸出是重建圖。 通過投影重建可以直接的看到原來被投影的物通過投影重建可以直接的看到原來被投影的物 體的某種特性的空間分布，比直觀觀測投影圖體的某種特性的空間分布，比直觀觀測投影圖 要直觀的多。要直觀的多。 5.6 圖像投影重建圖像投影重建 Radon變換 對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ù)和投影之間的關系，若函數(shù)為f (x, y)， 則不同角度下的投影可寫為 原理：原理：“斷層平面中某一點的密度值可看作這一平面內(nèi)所有經(jīng)過斷層平面中某一點的密度值可看作這一平面內(nèi)所有經(jīng)過 該點的射線投影之和（的平均值）該點的射線投影之和（的平均值）” 5.6 圖像投影重建圖像投影重建 Radon變換的matlab實現(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實現(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實現(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實現(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角度的投影值加進空白圖像 實 例 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 圖像投影重建圖像投影重建 反投影重建法反投影重建法 缺點：星狀偽影缺點：星狀偽影 000 010 000 原始圖像原始圖像 1/n1/n1/n 1/n11/n 1/n1/n1/n 重構(gòu)圖像重構(gòu)圖像 中心點中心點A經(jīng)經(jīng)n條投影線投影后，投條投影線投影后，投 影值均為影值均為1： p1=p2=.=pn=1 因此重建后因此重建后 而其他點均為而其他點均為1/n：這類偽跡成為：這類偽跡成為 星狀偽影星狀偽影 12 1 (.)1 An fppp n 5.6 圖像投影重建圖像投影重建 反投影重建法反投影重建法 缺點：星狀偽影缺點：星狀偽影 星狀偽影星狀偽影 5.6 圖像投影重建圖像

35、投影重建 反投影重建法反投影重建法 缺點：圖像模糊缺點：圖像模糊 圖像產(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、上沿同一方平面上沿同一方 向且過原點的直線上的值。向且過原點的直線上的值。 5.6 圖像投影重建圖像投影重建 傅里葉切片定理（中心切片定理）傅里葉切片定理（中心切片定理） 投影函數(shù)的數(shù)學表達式：投影函數(shù)的數(shù)學表達式： 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 , 目標函數(shù) f(x,y) 可由傅立葉函數(shù)F(u,v) 的逆 變換獲得，即 5.6 圖像投影重建圖像投影重建 濾波反投影法濾波反投影法 雅可比行列式 sin cos u 頻域中的笛卡爾坐標與 極坐標的關系為： 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(, )表示對應于角度的單位投影的傅立葉變換；里層的積分是 P(, )|的逆傅立葉變換，記為g(t,)，在空間域，它表示單位 投影被一頻域響應為|的函數(shù)做濾波運算，故稱之為濾波反投 影 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 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 頻域中：濾波器| 空間域中與其對應的濾波器為： det tj 2 將t0代入上式計算得到(0)，即曲線|以下的面積。 當時，(0)，所以上式是無法直接計算的，必須 另想它法，引入限帶函數(shù)（band-limiting function） 5.6 圖像投影重建圖像投影重建 濾波反投影法中濾波器的選擇濾波反投影法中濾波器的選擇 2jt ted 濾波器是個無

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

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

43、的選擇濾波反投影法中濾波器的選擇 Hamming濾波器 降低了高頻噪聲, 可得到Hamming濾波 器和Hanning濾波器 5.6 圖像投影重建圖像投影重建 反投影重建算法的反投影重建算法的matlab實現(xiàn)實現(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實現(xiàn)實現(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實現(xiàn)實現(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實現(xiàn)實現(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、況下的重建算 法較為復雜，但實質(zhì)法較為復雜，但實質(zhì) 沒有改變?？刹捎闷?jīng)]有改變?？刹捎闷?行束情況下的算法實行束情況下的算法實 現(xiàn)，只需加以適當?shù)噩F(xiàn)，只需加以適當?shù)?修正即可修正即可 n 重排算法重排算法: : 把一個視圖中采得的扇形數(shù)據(jù)重新組合成平行的把一個視圖中采得的扇形數(shù)據(jù)重新組合成平行的 射線投影數(shù)據(jù)，然后采用平行束重建算法重建射線投影數(shù)據(jù)，然后采用平行束重建算法重建 n 直接重建算法直接重建算法: : 不必數(shù)據(jù)重排，只需適當加權(quán)即可運用與平不必數(shù)據(jù)重排，只需適當加權(quán)即可運用與平 行束類似的算法重建行束類似的算法重建 5.6 圖像投影重建圖像投影重建 扇形束重建- EA等角扇形束重建

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

60、b實現(xiàn)實現(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