第二章 影像重組理論
2.4 濾波反投影演算法實作
FBP 演算法由前一小節可以知道，共分為濾波投影及反投影兩部分。根據 Filtered Projection Qθ
( )
t =∫
−∞∞Sθ( )
ω ωej2πωtdω (2.20) Back Projection f( )
x,y =∫
0πQθ(
xcosθ + ysinθ)
dθ (2.21)兩程式，假設 0°到 180°之間對待測物體共作K次不同角度的投影，所以演算法簡單的流 程如下：
I. 濾波投影
取得θi的投影資料Pθi
( )
t ，i從 0 到K −1，執行下列工作共K次。i. 將Pθ
( )
t 作傅立葉轉換到頻率域，成為Sθ( )
ω 。ii. 對Sθ
( )
ω 乘以濾波器ω在頻率域作加權。iii. 將加權後的Sθ
( )
ω 值作反傅立葉轉換回空間域(space domain)，成為Qθ( )
t 。 II. 背投影i. 從點
(
x,y)
及角度θi算出tj = xcosθi +ysinθi，i從 0 到K −1，共K次。ii. 累加每個Qθi
( )
t 中，位置tj上得到的Qθ( )
tji 值，共累加K次，得到
)
點原本的線性衰減係數。
(
x,y在實作上，電腦斷層掃描取得的投影資料會排列成一張影像，我們稱為sinogram。
sinogram的排列，是按照 0°到 180°依序將投影資料以一列一列的方式擺列成一張影像，
其中若平行投影光束有 256 束，則每列就有 256 筆偵測器得到的偵測資料。sinogram示 意圖如圖 (a)所示：
Ray beam width
0°
Projection Data
90°
180° (a)
(b) (c)
圖 2-7 (a)是 sinogram 影像的示意圖，(b)是 256x256 pixels 灰階的範例圖，(c)是 (b)的 sinogram
圖 2-7(b)是 256x256 的灰階影像，其中以座標(100,80)為中心點有個 3x3 pixels的小 正方形。而圖 (c)就是對圖 (b)進行 256 次投影，每次投影偵測到 256 筆偵測資料而得 到的sinogram影像，因此sinogram影像也是 256x256 pixels大小，在此我們重建影像要 求得的pixel點衰減係數就是pixel的灰階值。
我們參考前面的提到的 FBP 演算法流程，先在個人電腦上以 C Language 實作此演 算法，以驗證 FBP 演算法的正確性。根據演算法流程，FBP 演算法程式設計流程圖如下 所示：
Input
sinogram image Filter Backproject
Output reconstructed image
圖 2-8 FBP 演算法程式流程
所以程式上的內容主要也是分成四個部分，以底下的 pseudo code 程式來表示我們
軟體的實作內容：
double filtered_rel[angles][rays]; /* 儲存完成 filter 後的 singogra 影像 的陣列 */
double reconstructed[rays][ rays]; /* 儲存重建完成影像的陣列 */
main() begin
/******從 sinogram 影像中讀取投影資料******************/
read sinogram image to sinogram[angles][rays];
/******底下是filtered projection的pseudo code **************/
/* 作 1D FFT 及 IFFT 之前要作 2 冪次方的 ZERO Padding */
int logn = int(ceil(log10(rays)/log10(2))+1);
int fftsize = int(pow(2,logn));
/* 定義FFT及IFFT中存實部，虛部的array */
double *rel = new double [fftsize+1];
double *img = new double [fftsize+1];
/* 定義 window值 */
double *window = new double [fftsize+1];
for(i=1; i<fftsize/2; i++)
window[i+1] = window[fftsize-i+1] = i;
window[0] = window[1] = 1;
window[fftsize/2+1] = fftsize/2;
for( a=0 ; a<angle ; a++ ) loop
memset(rel,0,sizeof(double)*(fftsize+1));
memset(img,0,sizeof(double)*(fftsize+1));
memcpy(rel, sinogram[a],sizeof(double)*rays);
FFT(rel, img);
/* 做完 FFT 後與 windows 相乘 */
for(i=1;i<fftsize+1;i++) loop
rel[i]*=(window[i]/(double)fftsize*2);
img[i]*=(window[i]/(double)fftsize*2);
end loop;
IFFT(rel, img);
memcpy(filtered_rel[a],rel,sizeof(double)*rays);
end loop;
/******底下是back projection的pseudo code **************/
/* 建立三個tables 查表用 */
for(x=0; x<angles; x++) {
block_x_cen = block_size;
block_y_cen = block_size;
shift_x = sinValue;
shift_y = cosValue;
x_sino_oricoor = shift_mp - mp*shift_y ; y_sino_oricoor = shift_mp + mp*shift_x ; x_sino_endcoor = shift_mp + mp*shift_y ; y_sino_endcoor = shift_mp - mp*shift_x ;
u_vector_x = block_x_cen - x_sino_oricoor;
u_vector_y = block_y_cen - y_sino_oricoor;
v_vector_x = x_sino_endcoor - x_sino_oricoor;
v_vector_y = y_sino_endcoor - y_sino_oricoor;
first_block_sinopoint[x] = (int)(((u_vector_x * (v_vector_x/1024)) + (u_vector_y *
(v_vector_y/1024)))/(128));
mod_block_x_cen = block_x_cen ; mod_block_y_cen = block_y_cen ;
x_coor_array[x] = mod_block_x_cen + shift_out * shift_x + ((-1)*block_sino_shift*shift_y);
y_coor_array[x] = mod_block_y_cen + shift_out * shift_y - ((-1)*block_sino_shift*shift_x);
}
/*區塊式影像重建*/
for(block_up=0; block_up<x_max; block_up+=block_size) {
count_y++;
count_x=-1;
for(block_left=0; block_left<y_max; block_left+=block_size) {
block_right = (block_left +block_size);
block_down = (block_up +block_size);
count_x++;
{ ((shift_y*count_x - shift_x*count_y)*block_size);
block_sino_start = (block_sino_cen - block_sino_shift);
block_sino_end = (block_sino_cen + block_sino_shift);
x_coor = x_coor_array[x];
y_coor = y_coor_array[x];
for(delta_r＝block_sino_start to block_sino_end) {
if(out_x_coor & out_y_coor in the boundary) {
result[block_size*count_y+out_y_coor][block_size*cou nt_x+out_x_coor] += filtered_rel[x][delta_r];
}
write reconstructed[rays][ rays] to reconstructed imgae;
end;
上面的pseudo code內容是針對 128x128 大小的sinogram影像(0°到 180°共 128 個投 影角度，每次投影偵測有 128 筆偵測資料)可以重建出截面圖為 128x128 pixels大小的 的灰階影像。pseudo code中，在濾波投影(filtered projection)過程時，傅立葉轉換 我們以快速傅立葉轉換(Fast Fourier Transform，FFT)來完成，且執行快速傅立葉轉 換前，取樣資料需先執行Zero Padding擴大成 2 的冪次方。而在頻率域中需使用的濾波 器ω 函數，則如圖 所示，將投影資料的一維傅立葉值⎯Sθ
( )
ω ，乘與一個權重，低頻 中心點權重值低，越往高頻則權重越高。圖 2-9 頻率域中的濾波函數
第三章 區塊化重建影像
本章我們設計一個區塊式重建影像的概念，作為整章闡述的主體，經由實作證實去 我們設計的演算法，以下將逐步說明區塊化重建影像概念：
3.1 反投影影像重建概念
由於 ML310 所提供的硬體資源不足以保存整個重建結果，所以必須將重建過程分塊 處理，我們採用 4x4 的遮罩作為重建區塊，依序重建出結果，當每個區塊都重建完成後，
即為最後的重建影像。如圖 3-1 所示，原圖大小為 128x128，每個區塊大小為 4x4：
圖 3-1
利用 X-ray 穿透待測物之後得到的衰減值，對應的 x-ray 射線透射物體時直線經過 的點的衰減係數總和，如圖 3-2 所示，
圖 3-2 投影示意圖
相當於是 x-ray 射線經過物體路徑的線性積分。如圖 3-3 所示將P(θ0)~P(θ127)由 上而下排列即為 projection data（sinogram）:
圖 3-3
圖 3-4 為 sinogram，圖 3-5 所示為 filtered sinogram：
圖 3-4 sinogram 圖 3-5 filtered sinogram
圖 3-7 Backprojection 示意圖
3.3 演算過程推演
當光線直接反投影回去，直接將 累加在重建影像的像素中，此時
我們取整數點當作反投影時所經過的座標，當 向前傳播 時，X 方向
的變化量為
) , ( _real i j filtered
) , ( _real i j
filtered Δk
θ
×sin
Δk ，Ｙ方向的變化量為Δk×cosθ ，如此可以方便計算出 下一個將要經過的座標，
) , ( _real i j
filtered Δk每次遞增值為一，以影像為例，左上角為
原點，因此，如圖 3-8 所示：
圖 3-8
圖 3-9
我們採用區塊式做影像重建，所以每一個區塊只需截出一小段的 sub-sinogram 來 做反投影。我們先求區塊中心(即 a 點)、 (即 b 點)、
圖 3-10
3.4 拆解區塊式 FBP
我們預先將反投影所需的資訊編成訊息列（message list），每筆訊息長度為 64 位 元，訊息內容包括 、Ｙ座標、Ｘ座標，各佔 32 位元、16 位元、16 位 元。如圖 3-11 所示：
) , ( _real i j filtered
圖 3-11
將七個訊息分別編成上述格式之後，再以平行的方式掃過整個區塊，被掃過的區域 立即累加 值，如圖 3-12，由 A 端掃到 B 端，A 端與 B 端的距離為 8 個 單位長，每條訊息所走過的路徑長均為 8，如黑色訊息所掃過的區域共有五塊。
) , ( _real i j filtered
圖 3-12
建立三個表格存放重要的資訊，利用查表法可以更快速將重覆且複雜的計算簡化成 表格，用簡單的位移即可計算出其他區塊所需的資訊，我們用的圖大小為 128x128，每 個區塊大小為 4x4，因此每一列可以切成 32 個區塊，總共有 1024 個區塊，0≤ m≤31，
，我們建立的表格資訊皆由 block(0,0)計算得來的，有了 block(0,0)資訊，
即可以推算出 block(m,n)於重建時所需的資訊，可計算出 block(m,n)中心點對於 31
0≤ n≤
) ( i P θ 的投影點，進而截出反投影時所需的 sub-sinogam，不用每次重建某一個區塊時，就重 新計算 sub-sinogram，節省不少計算上的時間，如圖 3-13：
圖 3-13
所要建立的三個表格分別儲存αi(0,0)，0≤ i≤127，截出重建時所需要的 sinogram 之後，第一個點的 X 座標與 Y 座標，即圖 3-14 中標示紅色的點，αi(0,0)為 blcok(0,0) 在P(θi)上的投影點：
圖 3-14
每個區塊中心點對P(θi)的投影點關係，經由推導可得
圖 3-16
有鑑於上述情形的發生，我們將訊息列分成偶數訊息(even messages)與奇數訊息 (odd messages)，偶數訊息有自己的重建區塊，奇數訊息亦有自己的重建區塊，奇數訊 息和偶數訊息重建之後的結果點對點相加起來，即為最後的重建圖形，偶數訊息與奇數 訊息於反投影過程中不會有資料衝突的情形，如圖 3-17 所示：
圖 3-17
第四章 結果與討論
4.1 結果與討論
以下為軟體模擬的結果，利用 Xilinx 所開發的一套嵌入式系統開發軟體
EDK(Embedded Development Kit)進行軟體模擬，撰寫 C 語言，搭配 ML310 中的 powerPC 進行運算，時脈為 100MHz，圖 4-14、圖 4-15 分別為 128x128 和 256x256 的原始圖，圖 4-1~圖 4-6 影像大小為 128x128，圖 4-1 採用區塊大小為 4X4 的結果，圖 4-2 採用區塊 大小為 8X8 的結果，圖 4-3 採用區塊大小為 16X16 的結果，圖 4-4 採用區塊大小為 32X32 的結果，圖 4-5 採用區塊大小為 64X64 的結果，圖 4-6 採用區塊大小為 128X128 的結果，
明顯的圖 4-6 結果較好，因為不需要計算截出 sub-sinogram，可避免誤差累積，其他不 同大小的區塊重建結果，由於在截出 sub-sinogram 時，直接取其整數值，亦沒有做內 插作為些許的修正，因此結果中出現了假影(artifacts)，此演算法目標用於硬體加速 影像重建，我們所使用的硬體資源有限，最後我們採取區塊的大小為 4x4，此演算法可 高度平行化的做影像重建。
圖 4-1 圖 4-2 圖 4-3
圖 4-4 圖 4-5 圖 4-6
圖 4-7~圖 4-13 為影像大小 256x256 的軟體模擬結果,圖 4-7 採用區塊大小為 4X4 的結果，圖 4-8 採用區塊大小為 8X8 的結果，圖 4-9 採用區塊大小為 16X16 的結果，圖 4-10 採用區塊大小為 32X32 的結果，圖 4-11 採用區塊大小為 64X64 的結果，圖 4-12 採用區塊大小為 128X128 的結果，圖 4-13 採用區塊大小為 256X256 的結果，相同的道 理，圖 4-13 結果較佳，因此若有更多硬體資源的話，可得到較好的重建結果。
圖 4-7 圖 4-8
圖 4-9 圖 4-10
圖 4-11 圖 4-12
圖 4-13
圖 4-14 圖 4-15
4.2 未來與展望
本演算法可於 FPGA 硬體加速計算，可高度平行化處理計算，且可以利用 pipeline 排程，當每一列做完傅立葉轉換、乘上濾波器、反傅立葉轉換，即可立刻做反投影，平 行化處理加上 pipeline 排程，預期的效果可以更好，執行時間可縮短，更有效率，使 硬體效能發揮到極致。
參考文獻
[1] Alexandro M.S. Adario, Eduardo L. Roehe, Sergio Bampi,
“Dynamically Reconfigurable Architecture for Image Processor Applications＂, DAC99, New Orleans, Louisiana, ACM, 1999
[2] Didier LATTARD, Guy MAZARE, “Image reconstruction using an original asynchronous cellular array＂, ISCAS IEEE, 1989
[3] Anup B. Sharma, Keith R. Allen and Roy P. Pargas, “Some new systolic designs for two-dimensional convolution＂, ACM, 1988
[4] Bei-Chuan Chen, Yu-Tai Ching, “A new antialiased line drawing algorithm＂, Computers & Graphics, 2001
[5] Bertil Schmidt, Manfred Schimmler, Heiko Schroder, “Tomographic Image
Reconstruction on the Instruction Systolic Array”, Computers and Artificial Intelligence, 1995
[6] A.V. Lakshminarayanan, "Reconstruction from divergent ray data," tech. rep., Dept.Computer Science, State University of New York at Buffalo, 1975.