第二章 影像重組理論
2.1 影像重建演算法
2.1.3 濾波反投影演算法實作
FBP 演算法由前一小節可以知道，共分為濾波投影及備投影兩部分。根據 Filtered Projection f
( )
x,y =∫
0πQθ(
xcosθ +ysinθ)
dθ (2.20)Back Projection Qθ
( )
t =∫
−∞∞Sθ( )
ω ωej2πωtdω (2.21) 兩程式，假設 0°到 180°之間對待測物體共作 K 次不同角度的投影，所以演算法 簡單的流程如下：I. 濾波投影
取得θi的投影資料Pθ
( )
t ，i 從 0 到i K−1，執行下列工作共 K 次。
i. 將Pθ
( )
t 作傅立葉轉換到頻率域，成為Sθ( )
ω 。ii. 對Sθ
( )
ω 乘與濾波器ω 在頻率域作加權。iii. 將加權後的Sθ
( )
ω 值作反傅立葉轉換回空間域(space domain)，成為( )
t 。 Qθ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示意圖如圖 2-7(a)所示：
0°
90°
180°
Ray beam width Projection Data
(a)
(b) (c)
圖 2-7 (a)是 sinogram 影像的示意圖，(b)是 256x256 pixels 灰階的範例圖，(c)是(b)的 sinogram 影像。
圖 2-7(b)是 256x256 的灰階影像，其中以座標(100,80)為中心點有個 3x3 pixels
的小正方形。而圖 2-7(c)就是對圖 2-7(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 演算法程式流程
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 **************/
/* the mid-point是投影光束的中心位置 */
int midpoint = rays/2;
for(y=0; y<y_max; y++) loop for(x=0; x<x_max; x++) loop
double sum = 0;
for(int i=0; i<angles; i++) loop theta = i*PI/angles;
t =cos(theta)*((x-x_cen)/(double)(x_cen)) +sin(theta)*((y-y_cen)/(double)(y_cen));
/* mid_bin表示當角度為theta時，重建影像(x,y)pixel對應於Qθ(t)的位置*/
mid_bin = (t* mid_point)+ mid_point;
/* mid-bin可能在Qθ(t)的整數位置low-bin跟high-bin之間，所以要求
mid-bin在low-bin與high-bin之間的位置比例，方便以後作內插法 */
lb = (int)mid_bin;
hb = lb + (mid_bin>lb? 1:0);
frac = mid_bin - lb;
/* mid_bin可以能在low_bin跟high_bin中間，所以作gray level要作內插法 計算 */
sum += (1-frac)*filtered_rel[i][lb] + frac*filtered_rel[i][hb];
end loop;
reconstructed [y][x] = sum;
end loop;
end loop;
/******將重建好的截面圖pixel灰階值寫回重建影像中******************/
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 的冪次方。而在頻率域中 需 使 用 的 濾 波 器 ω 函 數 ， 則 如 圖 2-9 所 示 ， 將 投 影 資 料 的 一 維 傅 立 葉 值
⎯Sθ
( )
ω ，乘與一個權重，低頻中心點權重值低，越往高頻則權重越高。圖 2-9 頻率域中的濾波函數