由投影重建影象:濾波反投影、FDK、TFDK三維重建演算法理論基礎
[1] RafaelC.Gonzalez, RichardE.Woods, Gonzalez,等. 阮秋琦等譯.數字影象處理(第三版)[M]. 電子工業出版社, 2011.P232
[2] RafaelC.Gonzalez, RichardE.Woods, StevenL.Eddins. 阮秋琦譯.數字影象處理:MATLAB版:本科教學版[M]. 電子工業出版社, 2014. (第二版)P102
參考文獻[1]從P228-P245
瞭解:
1. 由投影重建影象
1. 計算機斷層(CT)的原理
2. 投影和Radon(雷登)變換(radontransform)
3. 正弦圖(sinogram)和Shepp-Logan幻影
4. 傅立葉切片定理
5. 使用平行射線束濾波反投影的重建
6. Ram-Lak濾波器與Hamming、Hann等窗函式
7. 使用扇形射線束濾波反投影的重建
1.由投影重建影象
利用Matlab來複現參考文獻[1]中的處理流程:
首先生成一幅原始影象:
%程式碼如下: I=zeros(512,512); for i=1:512 for j=1:512 if ((i-256)*(i-256)+(j-256)*(j-256))<1024 I(i,j)=1; end end end imshow(I,[]);
得到1°和90°兩個方向上的投影(一維),並回抹得到反投影影象(二維)
這兩個方向的反投影相疊加:
%程式碼如下:(這裡參考了MATLAB自帶的函式說明輸入:doc iradon開啟該參考說明) %就是怎麼得到特定角度下的一個投影,並形成反投影影象: R=radon(I,0:179); r1=R(:,2); II1=iradon([r1 r1],[1 1])/2; %得到1°角度下的投影影象 imshow(II1); r90=R(:,91); II90=iradon([r90 r90],[90 90])/2; %得到90°角度下的投影影象 imshow(II90); II1_90=II1+II90;%二者疊加 imshow(II1_90);
同理,得到45°和135°方向的反投影影象,並將這四幅反投影影象相疊加:
%程式碼如下:
r45=R(:,46);
II45=iradon([r45 r45],[45 45])/2;
imshow(II45);
r135=R(:,136);
II135=iradon([r135 r135],[135 135])/2;
imshow(II135);
II_1_90_45_135=II1+II45+II90+II135;%四個反投影影象相疊加!
imshow(II_1_90_45_135);
當反投影的角度取樣增多時,這裡從0°到179°每隔1°共採180個反投影,相疊加後形成反投影影象:(不使用濾波器時的反投影重建影象如下:)
%這裡用了iradon函式,線性插值,不使用濾波器。
ii2=iradon(R,0:179,'linear','none');
imshow(ii1,[]);
ii1=iradon(R,0:179,'linear','Ram-Lak');%預設使用線性插值,Ram-Lak濾波器!
該結果如下:
將重建後的影象與原影象進行對比:(左邊是原圖、右邊是重建得到的影象,注意模糊與振鈴現象)
我們再使用MATLAB自帶的一頭部幻影影象Shepp-Logan:
原圖:
未使用濾波器的重建影象:
使用Ram-Lak濾波器 :
使用Shepp-Logan濾波器:
%生成以上影象的程式碼如下:
P=phantom(512);
theta=0:179;
[R,xp]=radon(P,theta);
I1P=iradon(R,0:179,'linear','none');
I2P=iradon(R,0:179,'linear','Ram-Lak');
I3P=iradon(R,0:179,'linear','Hamming');
以下內容參考文獻[1]將濾波反投影的基礎知識過一遍:
笛卡爾座標系的一條直線由它的斜截式描述: ;或由其法線方程來描述:
平行射線束的投影可以由一組這樣的直線建模。如圖5.37所示,投影訊號中的任意一點由沿著直線 的射線和給出。連續變數的情況下,線求個變為線積分,由下式給出:
該式是沿xy平面內任意一條直線的f(x,y)的投影(線積分)的公式,就是雷登(Radon)變換。
符號R{f(x,y)}或R{f}有時用於代替(3)式中的 來表示f的雷登變換。雷登變換是由投影重建影象的基石,計算機斷層(CT)是其在影象處理領域的主要應用。在離散情況下,(3)式變為
正弦圖與雷登變換:
正弦圖包含了重建影象f(x,y)所需的資訊。
正弦圖的視覺分析僅限於實際應用、但有時對於演算法開發是有幫助的。
CT的關鍵目的是從投影得到物體的三維表示。其方法是反投影每一個投影,然後對反投影求和以產生一幅影象(切片),再堆積所有的結果以產生三維物體的再現。(這裡是指用二維切片堆積成三維體,與FDK/TFDK直接就是三維重建是不同的!)
5.11.4傅立葉切片定理
傅立葉切片定理即投影的一維傅立葉變換和被投影區域影象的二維傅立葉變換間的關係。
投影 的一維傅立葉變換為:
式(11)就是著名的傅立葉切片定理(或投影切片定理)。它說明了一個投影(一維)的傅立葉是得到這個投影的二維區域f(x,y)的二維傅立葉變換對應角度下的一個切片。正如圖5.41所示,任意一個投影的一維傅立葉變換可以沿著一個角度提取一條直線的F(u,v)的值來得到,而該角度就是投影時所用的角度。
下面推導濾波反投影公式,將用到傅立葉切片定理。
F(u,v)的反傅立葉變換為:
當c=0.54時,該函式稱為漢明窗(RichardHamming);
當c=0.5時,稱為韓窗(Juliusvon Hann)
加了窗函式的濾波器在空域的振鈴現象減弱。
我們可以預期,由於使用漢明窗的反投影有較小的振鈴,但稍微模糊一點。見下圖:
原圖 ;使用Ram-Lak濾波器 ; 使用Hamming窗加窗後的濾波器
在CT的多數應用中(特別是醫學上),像振鈴這樣的人為缺陷有嚴重的厲害關係,使其最小化是有意義的工作。調整濾波演算法(可以做文章的地方,一些碩士論文就這樣自創新的濾波器,發現效果有所改進,好,成文!比如文獻
[1] 張鑾. 基於平板探測器的錐束CT重建技術研究[D]. 中北大學, 2010.)、和硬體製造方面的改進(如,提供探測器的檢測細膩度,即取樣粒度)
因為斜坡濾波器(甚至在被加窗時)在頻率域的直流項為零,故每一幅反投影影象的均值將為零。這將意味著,每一幅反投影影象都將有正畫素和負畫素值,當所有的反投影影象相加形成最終的重建影象時,一些負畫素值位置可能變成正畫素,而平均值可能不為零,但是,典型地,最終的影象將還是有負畫素值。
當有關平均值的知識未知時,就使用標定的方法將影象的畫素值都歸一化到一個區間[0,255]。當典型的平均值的知識是可用的時候,可將該值加到頻率域的濾波器上,從而抵消斜波並防止直流項為零。當在空域中使用卷積時,截斷空間濾波器的長度(斜坡的反傅立葉變換)的真正效果都將防止其有零均值,這樣就完全避免的迫零問題。濾波反投影演算法:
1. 計算每一個投影的一維傅立葉變換;
2. 用濾波函式 乘以每一個傅立葉變換;
3. 得到每一個濾波後的變換的一維反傅立葉變換;
4. 對步驟3得到的所以一維反變換積分(求和)