影像重建演算法

第二章影像重組理論

2.1 影像重建演算法

影像重建演算法的主要問題在於如何利用投影資料來尋找物體截面圖上各個點的線性衰減係數。例如我們以x-ray對生物組織進行投影，則我們可以在投影面的另一邊偵測到x-ray的衰減量。在這裡我們考慮物體的每個點有不同的x-ray衰減係數(物質密度不同)，而每個偵測器得到的衰減量，就是對應的x-ray射線透射物體時直線經過的點的衰減係數總和(類比物質密度總和)，相當於是x-ray射線經過物體的那條路徑的線性積分(line integral)。如圖 2-2 所示：

圖 2-2

由上圖可知我們以二維函數 ^f

( )

^x^,^y ^{來描述物體，加上}

( )

θ^,^t 參數，我們以線性積分

使用 delta function

(2.2)

公式(2.3)就是函數的雷登轉換(Radon Transform)。集合所有同一個角度θ的值就成為角度θ的投影資料。來達到前述問題的封閉解(closed form solution)繼而重組出截面圖。連結傅立葉轉換 (Fourier Transform) 到待測物體的截面圖之間的基礎理論就是傅立葉切片理論 (Fourier Slice Theorem)。

2.2 Fourier Slice Theorem

連結傅立葉理論到投影資料與待測物體截面圖之間的基礎理論是由 Bracewell、

Ramachandran 和 Lakshminaraynana,發展，接下來的證明則是由 Kak 和 Slaney 推導的結果。二維傅立葉轉換定義如下：

(2.5)

過低頻中心點的射線。圖 2-3 可以清楚的表示公式(2.10)的意義

圖 2-3

有了上述結果，我們可以在θ¹，θ²，…，θ^k角度對物體作投影，得到投影資料後經過一維傅立葉轉換對應到截面圖的二維傅立葉頻率域中。如圖 2-4 所示：

圖 2-4

假設我們可以作無限個不同角度的投影，那我們得到的資料經過一維傅立葉轉換後，可以填滿截面圖的頻率域，接著我們只要再作反傅立葉轉換(inverse Fourier Transform) 就可以得到我們所希望的截面圖。

可是實際上，我們只能得到有限個不同投影角度的投影資料，且由圖可知，從投影資料經過一維傅立葉轉換對應到頻率域中的資料點是中心點(低頻)密集，外圍(高頻) 離中心點越遠越稀疏，因此頻率域沒有對應到的資料點就需要從已知資料點作內插法來當成預設資料。可是這其中會有相當大的誤差，尤其是越高頻的地方誤差越大，造成重建回來的截面圖會有影像剝蝕(image degradation)的情形。因為上面提到的是電腦斷層掃描重建影像概念上的模式，實作上我們得用不同的方式來完成。

2.3 濾波反投影演算法

這邊我們要使用來完成影像重建的演算法是濾波反投影 (Filtered Back Projection，FBP⁽¹⁾)演算法。這演算法也是由Fourier Slice Theorem衍生而來，是藉由公式原理的改寫來達成不同的運算實作。這邊我們先了解濾波反投影演算法的想法。

由前一小節我們可以知道角度θ的投影資料作一維傅立葉轉換可對應到物體截面圖二維傅立葉頻率域中的一條射線，這相當於截面圖二維傅立葉頻率域乘於一個簡單的濾波器(filter)而得到如下圖(b)所示的情形。

(a) (b) (c)

圖 2-5 這圖顯示從投影資料中得到的頻率域情形。(a)是理想狀況，(b)是經過傅立葉切片理論後得到的狀況，(c)是(b)圖經過濾波器加權後到的狀況，(c)圖中的資料狀況近似於(a)

但是對影像重建而言，我們希望得到的是投影資料的一維傅立葉轉換可以跟上圖(a)

的派外型射線是相當的，這樣我們集合幾個這樣的射線就可以重建出高品質的影像。對 Fourier Transform)定義如下：

(2.11)

將頻率域的標準座標系統轉成極座標系統(polar coordinate system)，如下 θ

(

ω^,θ +π

)

(

−ω^,θ

)

所以公式(2.21)稱為濾波投影(filtered projection)。而公式(2.20)則對每個Q_θ

( )

t 進行背投影(back projection)的動作。對於重建影像上的每個點

(

，在角度θ下可以找到對應的值

備註 1：後面的章節中，濾波反投影演算法將以 FBP 演算法來簡稱。

圖 2-6 背投影過程示意圖

2.4 濾波反投影演算法實作

FBP 演算法由前一小節可以知道，共分為濾波投影及反投影兩部分。根據 Filtered Projection ^Q_θ

( )

^t =

∫

−^∞∞^S_θ

( )

ω ω^e^j²^πω^t^dω (2.20) Back Projection ^f

( )

^x,^y ⁼

∫

0^π^Q_θ

(

^xcosθ ⁺ ^ysinθ

)

^dθ (2.21)

兩程式，假設 0^°到 180^°之間對待測物體共作K次不同角度的投影，所以演算法簡單的流程如下：

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

)

^及角度θⁱ^算出tj = xcosθi +ysinθi，i從 0 到K −1，共K次。

ii. 累加每個^Qθ_i

( )

^t 中，位置t^j上得到的Q_θ

( )

i ^{值，共累加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(θ₀)~P(θ₁₂₇)由上而下排列即為 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，可避免誤差累積，其他不

在文檔中區塊式濾波反投影演算法設計與FPGA實作(I) (頁 13-0)

第二章 影像重組理論

2.1 影像重建演算法

( )

( )

2.2 Fourier Slice Theorem

2.3 濾波反投影演算法

(

)

(

)

( )

(

2.4 濾波反投影演算法實作

( )

∫

( )

( )

∫

(

)

( )

( )

( )

( )

( )

( )

(

)

( )

( )

)

(

( )

第三章 區塊化重建影像

3.1 反投影影像重建概念

3.3 演算過程推演

3.4 拆解區塊式 FBP

第四章 結果與討論

4.1 結果與討論

第二章影像重組理論

第三章區塊化重建影像

第四章結果與討論