計算 MI 與影像對位的實作

volume data 記作 window，然後針對 target 圖中選取

target volume data，也就是 target voxel。圖 4.1 是 MI 整體流程的示意圖，

其中 MI 的計算是 MI(A,B)=H(A)-H(B)-H(A,B)，要計算 H(A)以及 H(B)要計算 1D histogram，要計算 H(A,B)則是要計算 joint histogram，所以我們在 CPU 及 GPU 的實作都是把 H(A)和 H(B)以及 H(A,B)分開計算，再求出最後的 MI 值。我們也 針對醫學影像的特殊 histogram 的分佈，做出了特殊的處理以利於 GPU 的平行運 算結構。

圖 4-1 MI 流程示意圖

我們用CUDA實作H(A)和H(B)包含了以下五個步驟:

1.在GPU上宣告A圖和B圖的記憶體空間，還有GPU上儲存entropy的記憶體空間 2.把A圖和B圖從CPU的記憶體上傳到GPU的記憶體空間

3.決定kernel的gridDim (要有多少個blocks)和blockDim(每個block有多少個thread) 4.執行kernel後entropy的結果會存在GPU的memory內

5.把entropy的結果從GPU的記憶體送回到CPU的記憶體 是 sequential 的去 search data 中找尋相同大小的 window，用函式的參數去控制空 間上的位置，循序的計算 MI 後再比較大小。在 CUDA 的實作中，我們是用一樣 的想法，但是利用每一個 block 去作一個 window，也就是 volume data，由於一 個 kernel 可以分成非常多個 block，我們也就可以平行的處理好幾塊 volume data 的 MI，也就是平行一次計算好幾塊的 MI，然後比較彼此間 MI 的大小，而 block index 可以控制空間上的位置。

由於桌上型電腦的GPU最少需要處理桌面的更新，如果顯示卡執行同一個 kernel時間過久，可能會使得 Windows 螢幕的畫面沒辦法更新。我們的做法是 變成一個 iterative algorithms，不會把所有的volume data都在同一個kernel內執 行，而是利用iterative的方式一次計算一些volume data，這樣也不必反覆的在GPU 和CPU之間傳輸影像資料。因為圖A和圖B是已經傳送到GPU裡面，可以一直給 kernel使用，我們只需要把entropy的結果從GPU傳回CPU就好。

在 kernel 的部份，這是 cuda 最重要的地方，kernel 也是實際在 GPU 上執行 的程式，我們想法是一個 grid 平行的處理很多個 target window，由一個 block 去 處理一個 targer window，我們分配給一個 block 32 個 thread，以利於我們利用 intra-warp 的方式去處理 shared memory 的 write collisions，並且以一個 block 去

處 理 一 個 window 我 們 就 能 用 __syncthreads 去 處 理 block level 的 syntronization，我們才能確定所有 thread 已經存好了 histogram，然後在 kernel 內把 entropy 算出來，我們的 kernel 大概是這樣規劃的:

H(A)和 H(B)的部份

1. 初始化 shared memory 空間以及設定好每個 thread 的 threadTag。

__syncthreads

2. 框出該 block 要處理的 volume data window，並且平行的由 block 中的 thread 去讀取圖的 intensity，利用 shared memory 儲存 histogram 的結果，藉由 4.2 敘述 的方式處理 shared memory 的 write collisions。至於 intensity 較集中的部份則是每 一個 thread 都有自己的 shared memory 空間去存取。

__syncthreads(同一個 block 裡面的 threads 都要完成了 histogram 的計算)

3. 將每個 thread 中儲存 intensity 是較集中的部份的 shared memory 空間分別累加 起來。

__syncthreads

4.利用計算出來的 histogram 去計算 Marginal Probability Distribution，進而利用 cuda 內建的_log 函式來求出 entropy value。

我們用CUDA實作H(A,B)的想法跟實作h(A)和H(B)有些類似，只是在處理 joint histogram 需要的記憶體空間非常大，例如說一個256 bins的影像如果算1d histogram只要用到256x4個bytes(大約1KB)去存，但是joint histogram就要用到 256x256x4個bytes(大約256KB)的記憶體空間，所以我們就要用到global memory 來儲存。global memory是有硬體支援的atomic函式的，由於shared memory的存入 跟讀取的時間較快，我們還是會將2D的joint histogram像素值比較集中的部分由 shared memory存入跟讀取，其他的部分再由global memory來存入。我們計算 H(A,B)分成兩個kernel去做，其中一個kernel把H(A,B)中像素值較集中的部份用 shared memory去存，另外一個kernel存取其他部份用global memory，由於joint

entropy的算法是H(A,B)= ，是分別計算各

去讀取 A B gray level pair shared memory histogram 4.2 敘述的方式處理 shared memory 的 write

collisions，至於 (0,0) 4.2 提及的一樣

thread 都有自己的 shared memory gray level pair (0,0)

__syncthreads(同一個 block 裡面的 threads 都要完成了 histogram 的計算)

3. 合併 gray level pair(0,0)的部份的 shared memory 空間。

__syncthreads

4.利用計算出來的 histogram 去計算 Probability density function，進而利用 cuda 內建的_log 函式來求出 entropy value。

計算 H(A,B)中的其他部分

1. 框出該 block 要處理的 volume data window，並且平行的由 block 中的 thread 去讀取 A 圖跟 B 圖中的資料，然後可以得到 gray level pair，如果 gray level pair 的像素值都不是 0 則利用 global memory 儲存 histogram 的結果，藉由 CUDA 內 建的 atomic 函式，可以計算出我們想要的 histogram。

。

__syncthreads(同一個 block 裡面的 threads 都要完成了 histogram 的計算)。

2.利用計算出來的 histogram 去計算 Probability density function，進而利用 cuda 內建的_log 函式來求出 entropy value。