中 華 大 學

中 華 大 學

碩 士 論 文

題目：使用 CUDA 實作立體影像渲染於醫學影像

Using CUDA implemented volume rendering method applying to medical images

系 所 別：生物資訊學系碩士班 學號姓名：M09620005 廖紹宇 指導教授：林志陽 博士

中華民國 九十九 年 八 月

摘 摘

摘 摘 要 要 要 要

電腦斷層掃描（CT）是現今醫學上常使用到的三維攝影技術，醫生能透過電腦 斷層影像來進行體內組織之診斷與評估，而立體渲染法（Volume rendering）是醫學 診療上常用的三維視覺化技術，透過立體渲染法將這些連續的二維資料以三維的方式 繪製在二維平面上，能提供相較於二維顯示更多的資訊。本論文將以立體渲染法實行 CT影像的三維重建，針對不同的CT值範圍來設計轉換函數，來重建出人體組織之真 實色彩，為立體渲染法添加光照特效，以及利用預先積分法改善立體渲染法在低取樣 頻率下的失真問題，並且利用已有的GPU加速環境，在計算統一設備架構（CUDA）

環境下進行編程，取得較真實且高品質的CT三維重建結果。

關鍵字 關鍵字 關鍵字

關鍵字:立體渲染法立體渲染法立體渲染法立體渲染法、、、、計算統一設備架構計算統一設備架構計算統一設備架構計算統一設備架構、、、、預先積分預先積分預先積分預先積分

ABSTRACT

Nowaday, computed tomography(CT) is a common medical imageing method for produce a stack of 2D images. 3D volume rendering of CT images is a very useful technique for assisting doctors’ disgnosis. This work propose a CUDA base volume rendering technique for 3D volume rendering. We designed the transfer function by the range of the CT-value, and given real human tissue’s colors. Light illumination can produce better rendering quality. We also implemented a pre-integrated volume rendering method for reducing the alias resulted by undersampling rate, and using CUDA to accelerate the speed for creating the pre-integrated’s lookup table.

Keywords: computed tomography, volume rendering, pre-integrated

致謝 致謝 致謝 致謝

首先，非常感謝指導教授 林志陽博士，在這碩士班的日子裡不斷的細心指導 我，帶領我進入程式語言以及影像處理的世界，讓我有能力完成這篇論文。感謝荊 宇泰老師與黃俊燕老師百忙之中抽空參加我的口試，並且給予我指導。

在碩士班的這段日子裡，謝謝家慈、立德、文斌學長姐與統文、國慶、育欣、

柏妮你們的陪伴與關懷，你們的關心讓我倍感溫暖。

特別感謝我的家人，在大學和研究所的這段時間雖然很少回家，但是，由於你 們的支持，我才能完成這段學業，謝謝你們。

廖紹宇 僅識於中華生資所 中華民國 99 年 8 月

目錄 目錄 目錄 目錄

摘 要 … … … i

A B S T R A C T … … … i i 致 謝 … … … i i i 目 錄 … … … i v 圖 目 錄 … … … v i 表 目 錄 … … … v i i i 第 一 章 緒 論 … … … 1

1 . 1 研 究 動 機 及 目 標 … … … 1

1 . 2 立 體 資 料 渲 染 法 … … … 1

1 . 3 轉 換 函 數 與 光 照 特 效 … … … 2

1 . 4 繪 圖 處 理 器 泛 用 運 算 … … … 3

第 二 章 三 維 影 像 顯 示 技 術 … … … 5

2 . 1 光 線 投 射 法 … … … 5

2 . 1 . 1 立 體 渲 染 法 積 分 式 … … … 5

2 . 1 . 2 光 線 投 射 法 … … … 6

2 . 2 C T 值 與 轉 換 函 數 設 計 … … … 9 2 . 3 預 先 積 分 立 體 渲 染 法 … … … 1 1 2 . 4 光 照 特 效 … … … 1 2 2 . 5 計 算 統 一 設 備 架 構 … … … 1 3 2 . 5 . 1 C U D A 硬 體 模 型 … … … 1 4 2 . 5 . 2 C U D A 程 式 模 型 … … … 1 5 第 三 章 方 法 … … … 1 8 3 . 1 光 線 投 射 法 的 實 作 … … … 1 8 3 . 2 轉 換 函 數 的 實 作 … … … 2 0

3 . 3 預 先 積 分 立 體 渲 染 法 的 實 作 … … … 2 1 3 . 4 光 照 特 效 的 實 作 … … … 2 3 第 四 章 實 驗 結 果 … … … 2 5 4 . 1 實 驗 環 境 … … … 2 5 第 五 章 結 果 與 未 來 發 展 … … … 3 6 參 考 文 獻 … … … 3 7

圖 圖 圖

圖目錄 目錄 目錄 目錄

圖 1 - 1 人 體 胸 腔 的 C T 影 像 … … … 2

圖 1 - 2 近 年 G P U 與 C P U 浮 點 數 運 算 效 能 比 較 … … … 4

圖 2 - 1 光 線 投 射 法 原 理 模 型 … … … 7

圖 2 - 2 取 樣 / 三 維 線 性 內 插 法 … … … 7 圖 2 - 3 各 個 組 織 器 官 的 C T 值 範 圍 … … … 1 0 圖 2 - 4 預 先 積 分 立 體 渲 染 法 … … … 1 2 圖 2 - 5 C U D A 硬 體 架 構 … … … 1 4 圖 2 - 6 C U D A 程 式 架 構 … … … 1 5 圖 3 - 1 光 線 投 射 法 程 式 流 程 圖 … … … 1 9 圖 3 - 2 依 B l o c k 大 小 分 割 影 像 … … … 2 0 圖 3 - 3 轉 換 函 數 的 使 用 者 介 面 … … … 2 1 圖 3 - 4 預 先 積 分 法 程 式 流 程 圖 … … … 2 2 圖 3 - 5 預 先 積 分 法 於 光 線 投 射 法 之 程 式 流 程 圖 … … … 2 3 圖 3 - 6 光 照 特 效 程 式 流 程 圖 … … … 2 4 圖 4 - 1 預 先 積 分 立 體 渲 染 法 結 果 … … … 2 7 圖 4 - 2 預 先 積 分 法 與 較 高 取 樣 頻 率 結 果 比 較 … … … 2 8 圖 4 - 3 立 體 渲 染 法 之 光 照 特 效 結 果 … … … 2 9 圖 4 - 4 全 身 的 轉 換 函 數 之 設 計 … … … 3 0 圖 4 - 5 骨 骼 的 轉 換 函 數 之 設 計 … … … 3 1 圖 4 - 6 肺 臟 的 轉 換 函 數 之 設 計 … … … 3 2 圖 4 - 7 腦 部 的 轉 換 函 數 之 設 計 … … … 3 3 圖 4 - 8 血 管 造 影 的 轉 換 函 數 之 設 計 … … … 3 4 圖 4 - 9 腹 腔 血 管 造 影 的 轉 換 函 數 之 設 計 … … … 3 5

表 表 表

表目錄 目錄 目錄 目錄

表 2 - 1 各 組 織 器 官 R G B 值 … … … 1 0 表 2 - 2 C U D A 各 記 憶 體 之 比 較 … … … 1 6 表 4 - 1 預 先 積 分 法 之 二 維 轉 換 函 數 計 算 時 間 … … … 2 6

第一章 第一章

第一章 第一章 緒論 緒論 緒論 緒論

1.1 1.1 1.1

1.1 研究動機及目標 研究動機及目標 研究動機及目標 研究動機及目標

現今的醫學診療系統中，若要進行非侵入式的人體內部組織檢查，必須透過醫 學攝影技術來進行，如電腦斷層掃描（Computed Tomography, CT）、核磁共振攝影

（Magnetic Resonance Imaging, MRI）等三維醫學攝影技術，這些顯影技術都可以將 人體內部拍攝成為一張張的影像，這樣的一組連續的影像也就構築了人體的三維資 料，醫生可以透過這些影像來進行體內組織之診斷與評估。而這些由三維醫學攝影 技術所得到的連續影像除了以傳統二維方式顯示外，還能以三維方式顯示，也就是 藉由立體資料渲染法（Volume Rendering）將連續的二維資料以三維的方式繪製在 二維平面上，這種三維顯示方式相較於二維顯示所能提供的資訊更多、更廣泛，因 此許多研究都致力於三維影像之顯示。

計算統一設備架構（Compute Unified Device Architecture, CUDA）為 nVIDIA 所 提出的 API，是有別於傳統 OpenGL 以及 DirectX，CUDA 讓使用者能夠透過 C 語 言來編程及操控顯示卡，並且透過顯示卡的硬體架構實現平行處理運算，以降低大 部分重複性的運算時間。本論文利用已有的 GPU 加速環境，基於 CUDA 的立體光 線投射法實作 CT 影像的三維重建，並且參考原始 CT 影像之 CT 值，給予各個人體 組織之真實色彩；另外，為 CUDA 環境下的立體投射法添加光照特效，以及利用預 先積分法改善立體光線投射法在低取樣頻率下的失真問題，取得較真實且高品質的 CT 三維重建結果。

1.2 1.2 1.2

1.2 立體資料渲染法 立體資料渲染法 立體資料渲染法 立體資料渲染法

立體資料（Volume Data）是一種三維的資料結構，常用來描述真實世界中的

物體之內外結構，醫學上典型的立體資料如電腦斷層掃描、核磁共振攝影，如圖 1-1 所示為人類胸腔的 CT 影像，CT 或 MRI 按照一定的比例如每毫米擷取一個切面資 料，再將這些連續的二維資料堆疊成三維的立體資料。

圖 1-1 人體胸腔的 CT 影像

立體資料於三維空間分割上最小的單位我們稱為體素（Voxel），在一般的立體 資料中，體素是一個純量值，它代表了物體中某一個組織的強度大小，如溫度、密 度或影像中的灰階值。立體渲染法（Volume Rendering）是將三維立體資料以二維 投影的方式轉為平面影像的方法，立體渲染法的方法又可分為間接立體渲染法

（Indirect volume rendering）以及直接立體渲染法（Direct volume rendering）；直接 立體渲染法即是以三維立體資料作成像，最為基礎及核心的演算法為光線追跡法

（Ray tracing）[3]及光線投射法（Ray casting）[1，2]，另外尚有足跡法（Splatting）

以及切扭因子分解法（Shear-warp factorization）。本研究所採用的光線投射法是直接 立體資料渲染法的方法之一，其原理是模擬 X-RAY 穿透物體衰減的物理現象，假

階值給予某一個顏色來呈現，這種一個強度值對應一個顏色值和透明度的過程，我 們稱之為轉換函數（Transfer function, TF）；CT 影像中，由於不同的人體組織對 X-RAY 的吸收程度的不同，因此不同的器官組織在 CT 影像中都有其灰階值範圍，

Buzug[9]以及余[10]等著作中，我們可以整理幾個組織器官的灰階值範圍，再從人體 解剖圖擷取各器官的真實的顏色，設計出對應的轉換函數。

轉換函數的使用是為了讓立體渲染法能提供使用者更多資訊作參考，但是根據 奈奎斯特信號取樣理論（Nyquist sampling theorem），在取樣頻率較低時，信號的取 樣會引起交疊（Aliasing）的現象，導致信號的失真，提高立體渲染法的取樣頻率雖 然可以有效解決此問題，卻也提高了渲染的運算時間，K. Engel[4]提出了預先積分 立體渲染法，透過對一維轉換函數的前處理，提高了渲染結果的品質。

在立體渲染法中加入光照特效能夠有效的提高渲染結果的真實度，而最常運用 在立體渲染法的光照模型為 Phong 著色法[11]，有別於原本的 Phong 著色法僅針對 光線在物體表面的行為，在立體渲染法中考慮到立體資料中每一個取樣點與週遭體 素的梯度關係，當光線行進在立體資料中，我們以取樣點為中心取得 X、Y 及 Z 軸 方向作梯度的計算，直接提供 Phong 著色法做計算。

1.4 1.4 1.4

1.4 繪圖處理器泛用運算 繪圖處理器泛用運算 繪圖處理器泛用運算 繪圖處理器泛用運算

近年來，3D顯示卡上的繪圖處理器（Graphics Processing Unit, GPU）運算效能 不斷的提升，就以浮點數的運算來說，GPU的運算效能就已經是CPU的數倍，如圖1-2 所示。而由於頂點著色器（Vertex shader）以及像素著色器（Pixel shader）的演進，

讓使用者可以藉由程式的撰寫來控制GPU的運算，因此有人利用來做非圖形運算的應 用，這種非圖形運算應用我們稱為繪圖處理器泛用運算（General-Purpose Computation on GPUs, GPGPU）。

圖1-2 近年GPU與CPU浮點數運算效能比較[5]

由於GPU的硬體架構設計，使得GPU具有SIMD（Single Instruction Multiple Data）

的平行運算能力，而nVIDIA所提出來的計算統一設備架構（Compute Unified Device Architecture, CUDA）[5]也是一種GPGPU的技術，使用者可以藉由CUDA來進行程式 開發，將一個問題拆解成多個相同、且彼此不相關的小問題，再把這些分送至GPU 的Stream processor進行運算。

第 第

第 第二 二 二章 二 章 章 章 三維影像顯示技術 三維影像顯示技術 三維影像顯示技術 三維影像顯示技術

光線投射法（Ray casting）是立體渲染法的方法之一，此方法相似於光線追跡 法（Ray tracing），光線追跡法的概念是：由螢幕上的每一個像素點，對螢幕內的物 體投射一條光線，計算每一條光線在物體上的落點，並求出光源對物體亮度的貢獻，

當光線碰到物體後會形成折射或反射，再對折射或反射的光線作追蹤，最後將各受 光點的亮度作累加取得最後的像素值，當螢幕上每一個像素點皆經過光線追跡的計 算過後，便可以在螢幕上繪出此物體的三維模型。然而光線投射法並不會考慮反射、

折射等較複雜的光學計算，也就是僅考慮光線直接穿越立體資料，累積光線經過立 體資料的取樣點，此種不考慮光學行為的方法稱為光線投射法。

2 22

2. .. .1 11 1 光線投射法 光線投射法 光線投射法 光線投射法 2

22

2.1.1 .1.1 .1.1 .1.1 立體渲染法 立體渲染法 立體渲染法 立體渲染法積分式 積分式 積分式 積分式

根據 Nelson Max[7]所提出的光學模型，在現實世界中，我們肉眼之所以能看 見物體，是因為陽光（或環境光源）照射物體後，物體吸收或發散光線後，進入我 們的肉眼，最後投射在視網膜上。而立體渲染法演算法便是以物體吸收與發散光線 的特性為基礎架構發展。這裡以單一條光線的行進來解說立體渲染法，當一條光線 投射穿越過物體，依著這條光線行進的路徑，我們可以從這條射線上取樣得到一組 離散的座標點。我們可以依據每一個取樣點的強度與射線行進的距離

D

，得到射線 最終所呈現的顏色，如公式（2-1）[6]所示:

( )

( )

( ) ^{( ( ( ) ) )}

∫ ^{− ∫ ′ ′}

= ^D c s x t e ^{s x t d t} dt C

t

0

0

r τ r

（2-1）

其中式子的右項

x r ( ) t

代表了穿越物體的光線，所發生的吸收與發散等行為，其中

c

為取樣點的強度，

τ ( ) s

為光線的衰減係數（透明度），我們稱此方程式為立體渲染

法積分式。另外，立體渲染法積分式可以利用黎曼總和將光線行進長度分割成 n 等

分（代表取樣點總數為

^{n =}  ^{D / d} 

），近似於立體渲染法積分式的方程式如公式

（2-2）[6]所示:

( )

∑ ∏

⁻

=

=

−

=

1

0 0

~

ⁱ

1

j

i n

i

i

A

C

C

（2-2）

上述結果便是立體渲染法積分式對於光線之吸收與發散行為之描述，2.1.2 章節為立 體渲染法中的光線投射法之演算法。

2 22

2.1.2 .1.2 .1.2 .1.2 光線投射法 光線投射法 光線投射法 光線投射法

我們可以將光線投射法的演算法分為四個步驟: 光線投射（Ray casting）、取樣

（Sampling）、著色（Shading）及混合（Compositing），下面我們將對每一個步驟做 說明。

1. 光線投射（Ray casting）：我們可以將光線投射以一個簡單的照相機（或眼睛）

模型呈現，如圖 2-1 所示，假設一條由照相機中心射出的光線經過螢幕上的 像素後，繼續前進穿越立體資料。在這個步驟上，我們會檢視經過像素點的 射線是否有穿越立體資料，如果射線並沒有與立體資料有所交接，我們就不 會對這條射線作後續的計算；如果這條射線穿越了立體資料，便紀錄下射線 與 Bounding Box 之最遠的交點，當我們對立體資料作取樣時，一旦取樣點超 過最遠的交點，便停止下一個取樣。

圖 2-1 光線投射法原理模型。光線由螢幕上的像素投射，與立體資料最近之交點 N1~N5，以及最遠之交點 F1~F5。

2. 取樣（Sampling）：射線行進在立體資料時，射線每前進一段距離我們就做一 次取樣，但射線並不一定完整的落在體素上，而可能穿越體素與體素之間，

如圖 2-2 所示，因此利用三維線性內插法（Tri-linear interpolation）來計算取 樣點的灰階值。

圖 2-2 取樣/三維線性內插法[2]

eye

image pixel viewing ray

voxel

sample point trilinear interpolation

3. 著色（Shading）：

轉換函數的設計，透過轉換函數我們可以定義 CT 影像中每一個灰階值對應 到一個色彩值和透明度。

光照特效，由轉換函數我們可以得到每一個取樣點的色彩，再經由模擬光線 照射物體的光學現象，計算得到在光線影響下之色彩。

預先積分法，為了減少失真的結果，考慮前後取樣點的關係，預先將各種不 同取樣點顏色配對組合，計算儲存為一張的二維轉換函數。

4. 混合（Compositing）：混合可以分為兩種方式，由後至前混合（Back-to-front, BTF）以及由前至後（Front-to-back, FTB）。BTF 起始於最遠之取樣點以及起 始色彩

C

₀，沿著視線方向對各個取樣點作色彩之混合，假設取樣點

i

的位置 從 n-1 至 0：

( )

^' ₁

'

=

_i

+ 1 −

_{i i+}

i

C A C

C

（2-3）

本研究中所使用的混合方式為 FTB，FTB 起始於最近之取樣點，沿著視線方 向作混合，取樣點

i

的位置從 0 至 n：

( )

(

_i

)

_i

i i

i i i i

A A A

A

C A C

C

' 1 '

1 '

' 1 '

1 '

1 1

−

−

−

−

− +

=

− +

=

（2-4）

'

C

i與

A

_i^'為已混合之新的色彩值與透明度，由取樣點

i

的色彩值

C

_{i 以及透明}

2 22

2.2 .2 .2 .2 CT CT CT CT 值與轉換函數 值與轉換函數 值與轉換函數 值與轉換函數設計 設計 設計 設計

電腦斷層掃描是利用單一軸面的 X-RAY 旋轉照射人體，由於不同的組織器官 對 X-RAY 的吸收程度不同，在利用電腦的三維技術重建斷層面影像後，各個組織 所呈現的亮度也不同；在 CT 影像中，我們稱這些影像的亮度值為 CT 值，其單位 是 Hounsfield unit（HU），在理論上，以水為基準 0 HU，CT 值的範圍是介於-1000 HU 至 1000 HU 之間，而空氣的 CT 值定為-1000 HU，但是實際的 CT 影像中，有些物 質的 CT 值可能會超過 1000 HU，CT 值的範圍會來到-1000 HU 至 3000 HU，大約 等同於 12 位元的存取大小，因此在儲存 CT 影像時，也就是每一個灰階值皆以 12 位元存取。

由於 CT 影像每一個體素的值都是介於-1024 HU 到 3071 HU 之間，因此我們使 用大小為 4096 的一維轉換函數來指定各灰階值所該呈現之顏色與透明度，而各個人 體組織器官的 CT 值範圍我們可以在 Buzug[9]以及余[10]等著作中整理獲得，如圖 2-3 所示，因為每一個組織器官的 CT 值都有所不同，因此我們可以依其定義的範圍 給予色彩，如骨骼的 CT 值約在 250 HU 以上，所以我們可以將 CT 值坐落在這範圍 的體素，給予的骨骼顏色。

另外，為了呈現真實的人體組織顏色，我們最好的選擇便是由現實的彩色人體 解剖圖中擷取各個器官之色彩值範圍，這裡我們以 Visible Human Project[8]做範本，

針對幾種常需要由 CT 攝影來觀察的器官做色彩擷取，如表 2-1 所示。

圖 2-3 各個組織器官的 CT 值範圍[10]

表 2-1 各組織器官 RGB 值

組織 R 值 G 值 B 值

骨骼（緻密骨）

2.3 2.3 2.3

2.3 預先 預先 預先 預先積分 積分 積分立體渲染法 積分 立體渲染法 立體渲染法 立體渲染法

由於轉換函數定義了每一個灰階值所對應到的色彩值與透明度，讓使用者可以

對不同的 CT 值範圍作色彩上的調整，藉此找出並且突顯使用者有興趣的區域。然 而當我們替立體渲染法載入一維的轉換函數時，其結果往往會出現失真的情況，主 要 的 原 因 來 至 於 取 樣 頻 率 ； 根 據 奈 奎 斯 特 信 號 取 樣 理 論 （ Nyquist sampling theorem），在取樣頻率較低時，信號的取樣會引起交疊的現象，而導致信號的失真，

為了減少失真對立體渲染法結果的影響，提高取樣頻率是可以提高最後影像的品 質，但是相對的也會提升立體渲染法的計算時間。

預先積分立體渲染法解決了低頻率取樣的問題，原先設計一維轉換函數時，我 們只考慮灰階值對應到的 RGBA 值，而預先積分立體渲染法則是先考慮了前後兩個 取樣點之間的影響程度，假設

S

f 為前一個取樣點的灰階值，

S

_{b為後一個取樣點的}

灰階值，從設計好的一維轉換函數中，我們可以找到兩個取樣點對應到的 RGBA 值，然後我們事先將

S

_f 到

S

_b的每一種配對組合關係以立體渲染法的積分式計算，

其結果儲存為一張二維轉換函數中，如圖 2-4，而 K. Engel[4]提供了近似方法，如 公式（2-5）

( ) ( ( ) ( ) )

( ) ( ( ) ( ) ) _ ^





 



 −

− −

−

=

− −

=

f b

f b b

f

f b

f b b

f

S T S S T

S d d

S S

S K S

S K S d d S S C

exp 1 ,

, , ,

α

^（2-5）

其中，

^K ( ) ^{S =} ∫

^S

^C ( ) ^{S dS}

0

，

^T ( ) ^{S =} ∫

^S

( ) ^{S dS}

0

τ

_，

_d

為兩個取樣點之間的距離，由

公式（2-5）我們可以得到不同

S

f 與

S

_b組合下的 RGBA 值，並且儲存在一張二維 的轉換函數表中，當我們在做取樣時，便可以依照前後兩個取樣點的灰階值，至二 維轉換函數中找出適當的 RGBA 值。

圖 2-4 預先積分立體渲染法

2.4 2.4 2.4

2.4 光照 光照 光照特效 光照 特效 特效 特效

在實行立體渲染法時，為了提高其結果的真實度，有時會藉由光線照明來讓渲 染結果更為自然，而最常被使用到的光照模型為 Phong 著色法，Phong 著色法將復 雜的光學行為簡化成環境光（Ambient light）、漫射光（Diffuse light）以及鏡射光

（Specular light），由這三種光線來設定物體對光線吸收及反射的程度，在考慮光源 與觀察者的位置關係後，以公式（2-6）計算:

( )

_{s s}

( )

ⁿ

d d a

a

I k I l n k I h n

k

I v v v v

⋅ +

⋅ +

=

（2-6）

其中，

I

_a、

I

_d、

I

_s各代表著光源的環境光、漫射光以及鏡射光的強度；

k

_a、

k

_d、

k

為物體對環境光、漫射光以及鏡射光反射常數，用來描述物體表面對光線反射的

兩個體素之差值，如公式（2-7）所示，相較於前兩個梯度算法，六鄰差值能得到較 正確，而且滑順的結果。

( )

( ) ( )

( ) ( )

( ) ( )

T

z y x G z

y x G

z y x G z y x G

z y x G z y x G z

G y G x z G

y x G

 







 







−

− +

−

− +

−

− +

 =





 



∂

∂

∂

∂

∂

= ∂

∇

1 , , 1

, ,

, 1 , ,

1 ,

, , 1 ,

, 1 ,

, ,

,

_（2-7）

另外，就物理學的角度，一個點的光照度會和該點到光源的距離平方成反比，這 就是光衰減的概念。因此，我們可以由公式（2-6）增添光衰減算式，如公式（2-8）:

2

* 1

cd bd I a

I = + +

（2-8）

其中

a

、

b

、

c

為光衰減係數，而

d

為點到光源的距離；若

a

等於1，

b

與

c

為0，

則會產生不隨距離衰減的光線；若

b

等於1，

a

與

c

為0，則產生隨距離線性衰減的光 線；若

c

等於1，

a

與

b

為0，則產生隨距離平方衰減的光線。

2.

2.

2.

2.5 55 5 計算統一設備架構 計算統一設備架構 計算統一設備架構 計算統一設備架構

近年來，隨著顯示卡的發展，GPU 的效能也不斷提升，應用範圍越來越廣泛，

這歸功於 GPU 架構上有利於處理密集型數據以及平行數據計算。在 GPU 裡的資料，

每一個元素執行相同一項指令，因此降低了複雜的流程控制；也因為能夠同時運行 大量資料以及擁有強大的運算能力，較能隱藏資料計算存取所造成的延遲問題。

過去要利用 GPU 做圖形應用或非圖形應用時，我們會使用圖形應用程式介面 如：OpenGL、Direct3D，程式設計師必須根據 OpenGL 或 Direct3D 提供的函式設計 演算法以解決問題，在程式編寫上以及記憶體管理上則缺乏彈性；為了解決這些問 題，nVIDIA 發表了新的硬體以及程式架構，稱為計算統一設備架構（Compute Unified Device Architecture, CUDA），藉由 CUDA 來分配管理 GPU 的資料平行計算；透過 C 語言與 CUDA 的函式庫來編寫，因為不用使用傳統圖形函式庫處理非圖形應用計

算，使得程式設計更為方便。

2.

2.

2.

2.5 55 5.1 .1 .1 .1 CUDA CUDA CUDA CUDA 硬體模型 硬體模型 硬體模型 硬體模型

CUDA 的硬體模型擁有一套 SIMD 的多處理器（Multi-processor, MP），如圖 2-5 所示，每一個 Multi-processor 裡包含了四種晶片記憶體：Registers、Shared memory、

Constant cache 以及 Texture cache。在每一個時程中，Multi-processor 裡的每一個串 流處理器（Stream processor, SP）同時執行相同指令。舉例來說，本實驗所使用的 NVIDIA Geforce 210GT 架構上擁有 2 個 Multi-processor,每一個 Multi-processor 包含 8 個 Stream processor，所以總共有 16 個 Stream processor。

2.

2.

2.

2.5 55 5.2 .2 .2 .2 CUDA CUDA CUDA CUDA 程式模型 程式模型 程式模型 程式模型

圖 2-6 表示 CUDA 的程式架構，CUDA 如同其他圖形應用程式介面 OpenGL

以及 Direct3D，允許程式設計師運用 C 語言來做程式的編譯。圖 2-6 裡的 Host 指的 是 CPU 端，Device 指的是 GPU 端，在 CUDA 的程式架構中，主程式是由 CPU 端 來執行，而遇到了資料平行處理的部份，就會將要在 GPU 端執行的程式編譯成 Device 能執行的程式，在交由 Device 執行。而這段 Device 執行的程式在 CUDA 裡 稱為核心（Kernel）。在 CUDA 中，GPU 是能同時執行大量執行緒（Thread）的計 算裝置，每一群執行緒被組織成一個區塊（Block），如圖 2-6 所示。Block 是 Thread 的集合，且經由共享記憶體同步地、有效地共享資料。Block 能夠進一步的組織為 網格（Grid）。不過 Thread 在不同的 Block 不能傳遞以及同步。

圖 2-6 CUDA 程式架構

表 2-2 CUDA 各記憶體之比較

記憶體 領域 存取限制 位置 速度 記憶體大小

暫存器（Register） 執行緒 讀/寫 On-chip 相當快 32KB/MP

區域記憶體

（Local memory） 執行緒 讀/寫 DRAM 慢 DRAM 大小

共享記憶體

（Shared memory） 區塊 讀/寫 On-chip 相當快 16KB/MP

全域記憶體

（Global memory） 網格 讀/寫 DRAM 慢 DRAM 大小

常數記憶體

（Constant memory） 網格 唯讀 DRAM

（cache） 相當快 64KB

紋理快取

（Texture cache） 網格 唯讀 DRAM

（cache） 快 DRAM 大小

假如資料經過分割後暫存在 Multi-processors，每一個 Multi-processor 會擁有一 個或多個區塊，一個區塊內又包含多個執行緒，有多少個執行緒以及區塊將執行指 令，則視使用者分配區塊以及網格的大小而定。有效地將資料分配在每一個執行緒 是很重要的，多個執行緒以及區塊有便於執行指令時隱藏記憶體存取延遲問題。

表 2-2 為 CUDA 裡各個記憶體之比較以及其特徵。以下是本實驗主要所用到的 幾個記憶體：

1. 暫存器（Register）：GPU 端的程式中，大多數的變數都是以暫存器的形式存 放，這些暫存器是 Block 中所有的 Thread 共享的，也就是如果每一個 Thread

視角改變時，就會重新載入矩陣到常數記憶體。

第三章 第三章

第三章 第三章 方法 方法 方法 方法

在這個章節中，將介紹我們所提出的三維顯示系統，包括了在 CUDA 環境下 的立體渲染法、轉換函數介面的使用、基於 CUDA 的預先積分法，以及光照特效。

3.1 3.1 3.1

3.1 光線投射 光線投射 光線投射 光線投射法 法 法的實作 法 的實作 的實作 的實作

圖 3-1 為在 CUDA 環境下的立體渲染法程式流程圖，由於 CUDA 能夠藉由 C

語言以及其延伸函式庫來操控管理 GPU，因此我們在 Host 端進行立體顯示的環境 設定，以及 GPU 的記憶體管理，在交由 Device 端做立體渲染法的處理。流程共分 為 Host 端與 Device 端，以下為 Host 端：

1. 載入影像檔，並以 3D Texture 的方式傳送至顯示卡中的 Texture memory。

2. 設定 Pixel buffer object，依視窗大小重新計算 Block 大小。

3. 建立轉換函數，並以 1D Texture 的方式傳送至顯示卡中的 Texture memory。

4. 利用 OpenGL 架設場景，處理照相機的位置矩陣，並將矩陣傳入顯示卡記憶 體中。當視角改變時，會再重新計算矩陣，重新寫入顯示卡記憶體中。

5. 開始執行光線投射法。

圖 3-1 光線投射法程式流程圖

最後輸出影像上的每一個像素皆為投射出的光線，我們能以 Block[16,16]的大 小將影像切割，如圖 3-2 所示，影像大小為 512²，以 Block[16,16]切割後便有 1024 個 Block，然後 CUDA 將依硬體上所能支援 Multi-processor 數量分批執行光線投射

Load volume data to 3D texture memory

Set block size

Load transfer function to 1D texture memory

Set camera position and transform matrix

Run ray casting

Ray set

Have intersection with data?

Sampling Yes

Is farthest point?

No

Write to PBO Yes

No

Draw Start

Transfer function

Front-to-back compositing

Host Device

法，計算出每一個像素的色彩，以下是 Device 端:

1. 根據旋轉矩陣計算照相機與立體資料的相對位置，並找出射線與 Bounding Box 的交點，如果射線沒有穿過 Bounding Box，則直接將背景顏色寫入 PBO 中。

2. 從射線與立體資料之最近交點開始由前至後取樣，而取樣點的數值會由 GPU 從取樣點相鄰之八個頂點去做三維線性內插，每一個內插出來的數值再根據 轉換函數，找到其相對應的顏色與透明度，在依公式（2-4）做混合，直到混 合最遠的取樣點之色彩值為止。

3. 混合後的數值直接存入 PBO，當所有的射線都寫進 PBO，直接輸出於螢幕。

圖 3-2 依 Block 大小分割影像

之間；黑色線條即為轉換函數曲線，當使用者點擊滑鼠左鍵移動游標，就可變更該 CT 值的透明度。另外，我們在 A 區裡提供了立體資料的直方圖統計，便於使用者 在調整轉換函數時參考。在透明度圖層的上方 B 區與 C 區為轉換函數之色彩值，B 區是介於-1024 HU 至 3071 HU 之間，而 C 區是介於 0 HU 至 400 HU，使用者可以 藉由左上角 F 區的調色盤選擇所需的顏色，在點擊滑鼠左鍵移動游標改變各 CT 值 的色彩值；在圖的右上方 G 區為原始 CT 影像的各個切面。E 區為我們自訂的轉換 函數，在此介面點擊滑鼠右鍵便能從選單中選擇我們自訂的轉換函數設定；另外，

選擇不同的轉換函數會改變 G 區原始 CT 影像的窗值表現，D 區為目前轉換函數之 窗值的範圍。

圖 3-3 轉換函數的使用者介面。

3.

3.

3.

3.3 33 3 預先 預先 預先 預先積分 積分 積分立體渲染法 積分 立體渲染法 立體渲染法的實作 立體渲染法 的實作 的實作 的實作

根據我們在章節 2.3 中提到的取樣頻率所造成的失真問題，我們藉由預先積分

的方式是先將一維的轉換函數，利用公式（2-5）轉換為二維的轉換函數；但是在 CT 影像中，灰階值的範圍在-1024 HU 至 3071 HU，也就是有 4096 個的灰階值對應 到不同的 RGBA 值，當一維轉換函數再經過計算得到二維轉換函數後，其大小也就

來到 4096 × 4096 個 RGBA 值。我們利用 CUDA 平行運算的方式，將預先積分法 的計算在 GPU 端完成，以下為 GPU 實作預先積分立體渲染法的步驟:

1. 在 GPU 中的 Global memory 取得一塊 4096 × 4096 × 4 的浮點數空間，作為儲 存二維轉換函數之用。

2. 載入一維轉換函數至 GPU 中的 Texture memory。

3. 將 Block 大小設定為 16 × 16，因此我們可以劃分出 256 × 256 個 Block。

4. 以公式（2-5）計算獲得新的二維轉換函數。

5. 將二維轉換函數保留在 Global memory，當我們在執行立體渲染法時，便以取 樣點以及下一個取樣點的灰階值組合，在二維轉換函數中找尋適當的 RGBA 值。

Init GPU memory

Load transfer function to GPU

Pre-integrated Set block size

Run ray casting

進一步時，便對位置

i

以及

i + 1

取樣，提供轉換函數參考，但射線每向前一步便就 做兩次的取樣會使得渲染時間的提高，因此，在取得色彩值及透明度後，便將取樣 點

i + 1

儲存到取樣點

i

，讓下一次取樣時，僅需要做

i + 1

的取樣。

圖 3-5 預先積分法於光線投射法之程式流程圖

3.

3.

3.

3.4 44 4 光照特效 光照特效 光照特效 光照特效的實作 的實作 的實作 的實作

一般在電腦圖學中，要利用 Phong 著色法實現光線照射的特效，需要取得物體

表面上每一個頂點的法向量，藉此描述光線在各頂點反射狀況；但是在立體渲染法 中，使用的並非是物體表面的法向量，而是立體資料中每一個灰階值或是取樣點相

Sample point i

2D transfer function

Is farthest point?

No

Yes Front-to-back compositing

i = i+1 Ray set

Draw Sample point i+1

鄰的點以公式（2-7）計算，所得到的梯度值用來作為向量讓 Phong 著色法參考並建 立其光學行為。梯度值的計算可以選擇在執行立體渲染法之前，計算好的梯度值事 先儲存在與原始影像相同陣列大小的三維紋理中，當我們執行取樣動作時，再從紋 理中以三維線性內插的方式取得，這個方式的缺點就是必須使用到更多顯示卡記憶 體空間在儲存梯度值，尤其醫學影像取像技術逐漸的成長，影像的解析也越來越高，

使用更多的記憶體空間在儲存梯度值會造成更多的不便。這裡我們簡化圖 3-1 的流 程圖，並且將光照特效加入流程圖中，如圖 3-5 所示，本研究中計算梯度值的步驟，

則是每一次立體渲染法在做取樣時進行，當光照特效啟動之後，每一次取樣我們便 以取樣點為中心對各軸，依原始影像大小比例取得相臨六個點之灰階值，經公式

（2-7）計算得到取樣點之梯度值。

Start

Set light source and transform matrix

Run ray casting

Ray set

Sampling and Transfer function

Sample point’s gradient and relationship with light source

Front-to-back compositing

第四章 第四章 第四章

第四章 實驗結果 實驗結果 實驗結果 實驗結果

4 44

4.1 .1 .1 .1 實驗環境 實驗環境 實驗環境 實驗環境

以下是本實驗之軟硬體設備：

1. CPU：Intel Core2 Due 2.99GHz 2. 記憶體：3.25GB

3. GPU：nVIDIA GeForce 210

顯示卡記憶體：512MB

Multi-processor 數量：2 （共 16 Processor）

4. 作業系統：Microsoft Windows XP Professional 5. 開發軟體：Microsoft Visual .net 2005

6. 繪圖函式庫：OpenGL

GLUT：藉由 GLUT 來建立影像輸出視窗，利用滑鼠控制視角，另外開 啟轉換函數之視窗，以滑鼠及鍵盤調整顏色等事件。

Pixel Buffer Object（PBO）：計算過後的光線投射法結果寫進 PBO，在由 顯示卡將 PBO 輸出到螢幕上。

在本研究中，以 512²作為輸出影像大小之初始設定，渲染一張影像時間約為 91 ms~125 ms，每秒更新畫面數（Frames per second, FPS）大約為 8.5 fps~10.0 fps，渲 染時間的差異是因為不同視線角度下，實際穿越立體資料的射線數量，以及每一條 射線行進距離所引起的。

在表 4-1 中，我們可以看到以 CPU 計算預先積分法，並且將轉換函數載入顯示 卡中的 Texture memory，平均計算以及載入的時間約為 1172 ms，假設我們頻繁的更 新轉換函數，例如以手動方式在我們提供的介面進行顏色或透明度的微調，其渲染 結果將有所停頓才能呈現；我們使用 CUDA 執行預先積分法，取得大小 4096 × 4096

的二維轉換函數，其平均計算以及載入時間約為 109 ms。

表 4-1 預先積分法之二維轉換函數計算時間

平均計算時間以及載入轉換函數（ms）

CPU 端 1172 ms

GPU 端 109 ms

圖 4-1 為預先積分法之結果，其中 A、C、E 圖分別為未開啟預先積分之頭部

CT、胸腔 CT 以及胸腔 CT 局部特寫，B、D、F 圖則為開啟預先積分法之效果。在 A、C、E 圖中可以看見低取樣頻率下渲染結果有嚴重的水波紋影響，而 B、D、F 圖中的水波紋經由預先積分法的使用，大幅降低了水波紋對渲染結果之影響。

圖 4-1 預先積分立體渲染法結果。圖 A、C、E 為未開啟預先積分法之渲染結果，

圖 B、D、F 為開啟預先積分法之渲染結果。

在圖 4-2 中為我們比較在較低的取樣頻率使用預先積分法，與在較高的取樣頻 率不使用預先積分法之結果，其中 A 與 C 的取樣距離為每 0.01 做一次取樣，而 B 與 D 的取樣距離為每 0.005 做一次取樣，可以看見較低取樣頻率下使用預先積分法 的渲染結果與高取樣頻率渲染結果相近。

圖 4-2 預先積分法與較高取樣頻率結果比較。A 與 C 為低取樣頻率下開啟預先積分 法之渲染結果，B 與 D 為高取樣頻率下之渲染結果。

圖 4-3 立體渲染法之光照效果結果。A 為未開啟光照特效之原圖，BCDE 為不同光 源位置，分別為左方、後方、右方以及前方。

在章節 2.2 中，根據所獲得的各器官 CT 值範圍以及色彩值，我們另外提供幾 種不同情況之轉換函數設計:

全身（Body）:為我們設計的初始轉換函數，其範圍從-200 HU 至 2000 HU，低 於-200 HU 大多為空氣或其他低密度物質，高於 2000 HU 除了少數骨骼外，還有如 金屬等高密度物質，我們以脂肪、軟體組織、海綿骨、緻密骨等範圍作為主要顏色 設計。

予骨骼之 RGB 值，如圖 4-5 所示。

圖 4-5 骨骼的轉換函數之設計。

肺部（Lung）:肺部 CT 經常用來觀察肺臟中空氣的變化，而正常充氣的肺臟，

其 CT 值範圍在-900 HU 至-500 HU 之間，由於在這個觀察範圍內包含了空氣，因此 轉換函數的設計上，越接近 CT 值-900 HU 就以藍色作為空氣的顯示，越接近-500 HU 的 CT 值就以脂肪作為肺臟的顏色，如圖 4-6 所示。

圖 4-6 肺臟的轉換函數之設計。

腦部（Brain）:大腦白質的 CT 值大約是 20 HU 至 35 HU，大腦灰質得 CT 值大 約是 35 HU 至 45 HU。依照表 2-1 給予骨骼之 RGB 值，如圖 4-7 所示。

圖 4-7 腦部的轉換函數之設計。

CT 血管造影（Angiography）:血管造影是透過顯影劑來突顯血管與其他組織的 CT 值，幫助確定血管結構和術前規劃，而顯影劑的使用使得血管的 CT 值高於各器 官，這裡我們以兩種不同 CT 範圍來設計轉換函數，一個是以 100 HU 至 250 HU 作 為顯影劑的強度，如圖 4-8 所示，CT 值超過大部份的軟體組織，因此較適合觀察血 管；另一個範圍以 45 HU 至 250 HU，如圖 4-9 所示，腎臟、脾臟的 CT 值在顯影劑 的影響下，CT 值大約會來到 100 HU 以上，但是肝臟對顯影劑吸收的程度不如腎臟、

脾臟那麼強烈， CT 值僅提升約 15~20 HU，因此我們讓 45 HU 至 100 HU 以肝臟 顏色呈現。

圖 4-8 血管造影的轉換函數之設計。

圖 4-9 腹腔血管造影的轉換函數之設計。

第五章 第五章

第五章 第五章 結論與未來 結論與未來 結論與未來 結論與未來發展 發展 發展 發展

在本研究中，我們在CUDA環境下利用立體渲染法建立了CT影像的三維顯示系

統，為了增進立體渲染法的品質，我們結合了預先積分法，在較低的取樣頻率下也能 得到較好的渲染結果，並且透過CUDA的運算，使我們能較快的取得預先積分法的二 維轉換函數，使用光照特效，讓渲染結果更加的立體且真實。

在轉換函數的設計中，會遇到不同組織器官的CT值範圍相同，這會使得渲染結 果不易分辨出不同的器官，未來將加入影像的分割以及辨識來改善這種情況。

由於CT影像的取像技術越來越好，使得影像的資料量也跟著提升，有些CT影像 的資料量已經超過顯示卡的Texture memory的容量了，未來我們將設計適當的演算法 來處理這種較大的CT影像。在立體渲染法中加入光照特效使得結果更為真實，但在 本研究中並沒有考慮光線照射物體後產生陰影的結果，未來我們將在CUDA環境下針 對陰影的應用設計演算法，在不影響立體渲染法運行的情況下，有效的呈現完整的照 明效果。

參考文獻 參考文獻 參考文獻 參考文獻

[1] Roth, Scott D., "Ray Casting for Modeling Solids," Computer Graphics and Image Processing, No. 18: pp.109–144, 1982.

[2] Marc Levoy, "Display of surface from Volume Data " IEEE Computer Graphics and Applications, 1988.

[3] Glassner A. S., "An Overview of Ray Tracing," An Introduction to Ray Tracing, pp.1-31, 1989.

[4] K. Engel, M. Kraus, T. Ertl, "High-quality pre-integrated volume rendering using hardware-accelerated pixel shading," SIGGRAPH/Eurographics Workshop on Graphics Hardware 2001, pp. 9-16, 2001.

[5] Compute Unified Device Architecture (CUDA) Programming Guide.

[6] M. Hadwiger, J. Kniss, K.Engel and C. Salama, "High-Quality Volume Graphics on Consumer PC Hardware," ACM SIGGRAPH course, 2002.

[7] Nelson Max, "Optical models for direct volume rendering," IEEE Transactions on Visualization and Computer Graphics, pp.99–108, June 1995.

[8] MJ Ackerman, "The visible human project," Proceedings of the IEEE, 1998.

[9] T. Buzug, "Practical Aspects of Computed Tomography," Springer Berlin Heidelberg, pp. 471-484.

[10] 余建明, """"醫學影像技術學醫學影像技術學醫學影像技術學,"醫學影像技術學"""合記圖書出版,2009 年 7 月初版.

[11] B. T. Phong, "Illumination for computer generated pictures,"

Communications of the ACM, pp.311-317, 1975.