國
立
交
通
大
學
多媒體工程研究所
碩
士
論
文
磁核共振影像之腎功能
The study of kidney function in MRI scans
研 究 生：黃怡敏
指導教授：荊宇泰 博士
磁核共振影像之腎功能
The study of kidney function in MRI scans
研 究 生：黃怡敏 Student：Yi-Min Huang
指導教授：荊宇泰 Advisor：Yu-Tai Ching
國 立 交 通 大 學
多 媒 體 工 程 研 究 所
碩 士 論 文
A ThesisSubmitted to Institute of MultimediaEngineering College of Computer Science
National Chiao Tung University in partial Fulfillment of the Requirements
for the Degree of Master
in
Computer Science
October 2007
Hsinchu, Taiwan, Republic of China
磁核共振影像之腎功能
學生：黃怡敏
指導教授：荊宇泰 博士
國立交通大學多媒體工程所
摘
要
磁核共振造影可以提供影像比較好的空間解析度與時間解析度，是目前計算腎功能 的新方法，而計算腎功能時需要皮質的遮罩以及皮質範圍的灰階值變化，因拍攝時腎臟 會隨著呼吸而上下移動，所以我們將分別利用兩種方法將兩兩相鄰的影像的皮質部份做 校正，之後比較這兩種方法它們之間的優缺點，然後利用半自動的方式找到皮質的部 份，最後利用電腦協助計算出腎功能；此外，我們再將計算出來的腎功能與放射科醫生 所計算的值做比較，接著再與利用核子醫學所算出來的腎功能做比較，並探討它們之間 什麼情況之下會有不同的結果。
The study of kidney function in MRI scans
Student:
Yi-Min
Huang
Advisor:
Dr.
Yu-Tai
Ching
Institute of Mutimedia Engineering
National Chiao Tung University
ABSTRACT
Magnetic resonance imaging (MRI) is a new approach at present for calculation of kidney function by providing better Spatial Resolution and Temporal Resolution. The calculation of kidney function needs the cortex mask and the change of grayscale value in cortex. Due to the reciprocation of kidney with respiration, we use two approaches
respectively to align the region of cortex in each pair of contiguous images and compare the advantage and drawback between the two approaches. After that we use a semi-automatic method to find out the cortex. Then, calculate the kidney function via computer. Furthermore, we compare the kidney functions calculated by our approach with those calculated by
Radiologist. Finally, compare them with those calculated by Nuclear Medicine and discuss the situation causing different results.
誌 謝
本論文能得以完成，首先要感謝我的指導教授荊宇泰教授，在這兩年多來由於他的 專業知識傳授與細心的指導，才能讓我順利完成此論文。同時也感謝口試委員張北葉醫 師與陳永昇老師指正論文的謬誤處及提供許多專業上的建議，使本論文能更加正確完 整。感謝醫學影像實驗室的全體成員，在這兩年給予許多專業知識上及精神上的支持。 最後要感謝我的家人，感謝他們在這兩年多來給我的包容與支持，讓我得以安心的 完成我的碩士論文，謝謝大家。目 錄
中文提要 ……… i 目錄 ……… ii 表目錄 ……… iv 圖目錄 ……… v 一、 緒論 ……… 1 1.1 簡介 ……… 1 1.2 研究動機 ……… 2 1.3 論文架構 2 二、 背景與理論基礎 ……… 3 2.1 醫學背景……… 3 2.1.1 泌尿系統之代謝……… 3 2.1.2 腎功能之計算……… 4 2.2 傅立葉轉換 ……… 5 2.2.1 離散傅立葉轉換……… 5 2.2.2 快速傅立葉轉換……… 5 2.3 高斯低通濾波器 ……… 7 2.4 Canny 邊緣偵測……… 7 2.5 影像匹配 ……… 9 2.6 次影像 ……… 11 三、 實驗流程 ……… 13 3.1 第一種演算法 ……… 14 3.1.1 流程圖……… 14 3.1.2 針對時間軸做快速傅立葉轉換……… 15 3.1.3 利用高斯濾波器讓資料變化平滑……… 173.2 第二種演算法 ……… 23 3.2.1 流程圖……… 23 3.2.2 建構遮罩……… 24 3.2.3 資料前處理……… 25 3.3 兩種演算法之比較 ……… 27 3.3.1 MR urographic 之比較……… 27 3.3.2 時間之比較……… 32 四、 實作與結果 ……… 30 4.1 實驗環境 ……… 30 4.2 MRU 電腦協助計算與放射科醫師結果比較……… 32 4.2.1 放射科醫師之作法……… 34 4.2.2 比較……… 35 4.2.3 討論……… 36 4.3 MRU 與核子醫學結果之比較……… 37 五、 結論與未來展望 ……… 43 參考文獻 ……… 44
表目錄
表 3-1 實驗資料 ……… 13 表 3-2 資料一兩腎比例之比較 ……… 28 表 3-3 資料二兩腎比例之比較 ……… 30 表 3-4 資料三兩腎比例之比較 ……… 31 表 3-5 資料四兩腎比例之比較 ……… 32 表 3-6 資料五兩腎比例之比較 ……… 33 表 3-7 兩演算法時間之比較 ……… 34 表 4-1 誤差比例之分配表 ……… 36圖目錄
圖 2-1 MR urography ……… 3 圖 2-2 腎解剖示意圖 ……… 3 圖 2-3 腎功能計算示意圖 ……… 5 圖 2-4 高斯低通濾波器圖 ……… 7 圖 2-5 次像素計算示意圖 1 ……… 12 圖 2-6 次像素計算示意圖 2 ……… 12 圖 2-7 次像素計算示意圖 3 ……… 12 圖 3-1 實驗資料圖 ……… 14 圖 3-2 演算法一流程圖 ……… 15 圖 3-3 原始資料示意圖 ……… 16 圖 3-4 經高斯低通濾波器 ……… 17 圖 3-5 高斯參數 10 反傅立葉之結果 ……… 18 圖 3-6 高斯參數 20 反傅立葉之結果 ……… 19 圖 3-7 高斯參數 30 反傅立葉之結果 ……… 20 圖 3-8 高斯參數 40 反傅立葉之結果 ……… 21 圖 3-9 MRU 曲線之結果 ……… 23 圖 3-10 演算法二流程圖 ……… 24 圖 3-11 左右腎之遮罩圖 ……… 25 圖 3-12 前處理後 ……… 26 圖 3-13 有無前處理及次像素之比較圖 ……… 27 圖 3-14 資料一之右腎比較 ……… 28 圖 3-15 資料一之左腎比較 ……… 28 圖 3-16 資料二之右腎比較 ……… 29 圖 3-17 資料二之左腎比較 ……… 29 圖 3-18 資料三之右腎比較 ……… 30 圖 3-19 資料三之左腎比較 ……… 30圖 3-20 資料四之右腎比較 ……… 31 圖 3-21 資料四之左腎比較 ……… 31 圖 3-22 資料五之右腎比較 ……… 32 圖 3-23 資料五之左腎比較 ……… 32 圖 4-1 放射科醫師圈選之皮質 ……… 35 圖 4-2 放射科醫師使用之 MRU ……… 35 圖 4-3 電腦與放射科醫師結果之比較 ……… 36 圖 4-4 核醫錯誤之範例 ……… 37 圖 4-5 MRU 與核醫結果之比較 ……… 37 圖 4-6 MRU 與核醫相差值與腎長度之相關圖 ……… 39 圖 4-7 長度差(%)低於 40 之圖 ……… 40 圖 4-8 長度差(%)高於 40 之圖 40 圖 4-9 核醫偵測不到腎功能之範例 41 圖 4-10 MRU 與核醫相差值與皮質比例之相關圖 42
一、緒 論
1.1
簡介
腎臟是屬於泌尿系統的一部分，它的功能主要是過濾並且代謝有害物質， 亦可以維持體內電解質之平衡，一旦腎臟的腎功能發生問題，便會嚴重影響我們 的生活品質，若我們能越早知道自己腎功能的情況便能及早治療並挽回健康，因 此判斷腎臟的功能是否正常的確是不容忽視的問題。隨著醫學技術及醫學儀器的 進步，非侵入式的診斷造影(diagnostic imaging)安全又快速即成為協助醫師們診 斷的方法，而核磁共振與磁振造影(Magnetic Resonance Imaging, MRI)、電腦斷層 掃描(Computer Tomography, CT)、核子醫學造影等方法皆被廣泛使用於臨床用途。其中患有先天性泌尿系統擴張的病人需要功能及解剖資料的証明，傳統上
提供證明的方法為超音波造影(Ultrasound, US)、核素利尿性腎顯像(Diuretic Renal Scintigraphy, DRS)、泌尿圖(urography)，及排尿中膀胱尿道造影(Voiding
Cystoureterography, VCUG)，而 MRI 則可以同時為腎臟及泌尿管造影，於是便成 為目前最常使用的造影技術。 有了好的造影技術之後，我們希望能將擷取的影像資訊利用影像處理的方式 將其轉換為數據；而此論文內容主要改善在造影過程中，因呼吸而上下移動的影 像，我們利用兩種方法分別去校正它，比較他們之間的優缺點，接著利用影像切 割的方式去擷取我們感興趣的部份，最後根據醫學上的公式自動化計算腎功能， 最後統計出關於核子醫學所得之腎功能與 MRI 所求得之腎功能之間的差，在腎 水腫嚴重程度上的相關度。
1.2
研究動機
目前以核子醫學造影最常被醫師們當作診斷腎功能的依據，主要原因在於 計算腎功能時幾乎不需要人為的介入，故有一定的可信度，且價格較低；但自從 有了 MRI 技術之後，發現核子醫學造影可能在某些情況之下會發生計算錯誤， 例如左腎之腎功能從影像判斷後明顯高於右腎，但使用核子醫學造影計算後，結 果卻不如預期，因此我們希望藉由結合動靜態(combined static-dynamic)的 MRI 技術來做計算，與核子醫學影像造影之結果做相關比較，並分析可能發生錯誤的 情況。
1.3
論文架構
本篇論文共分為五個章節：第一章為緒論，簡介醫學影像的種類，及研究 動機與目的；第二章為背景與理論基礎，先會對得到的影像資料稍作說明，接 著將介紹本論文所使用到的基礎理論；第三章為實驗流程，共有兩個演算法， 第一個是針對在頻域領域的處理方法，第二個方法則是利用匹配的方式作校正 影像，最後為兩種演算法的結果做比較；第四章為實驗結果，因為第一種演算 法效率較高，所以此處則是利用第一種演算法處理過後，利用大量病歷資料算 出其數據，再利用數據做出統計表並且分析其分布情況，希望能夠協助臨床的 診斷；最後一章則為結論與未來研究的方向。二、背景與理論基礎
2.1
醫學背景
在計算腎功能時，我們需要有醫學的根據才能夠去計算其中所需要的參數及 資料，以下是關於計算腎功能所需要的背景。 2.1.1泌尿系統之代謝
在一組動態 MR urographic 資料中，如圖 2-1，它的影像平均灰階值變化如 下：在人體注射對比劑一段時間後，因為對比劑會隨著血液到達腎小球，使得 在腎臟周圍皮質(cortex)部分灰階值逐漸增大，此時反應出皮質組織從血液攝取 對比劑的能力，接著對比劑會從皮質代謝至髓質(medulla)，再到腎盂(pelvis)， 最後至腎臟外的收集系統(collecting system)；由此我們可以從對比劑的代謝過 程得到收集系統和泌尿系統的型態學影像(morphological image)。 圖 2-1 MR Urographic圖 2-2 腎解剖視意圖
2.1.2
腎功能之計算
我們根據[4]的方法來計算各別腎功能(split renal function)，方法主要是基 於在腎臟我們感興趣的部份(regions of interest, ROIs)所形成的時間-灰階值的 曲線圖。 然而在正常情況的腎臟中，MRU 的腎圖中可以觀察出有三種典型的階 段：第一階段，曲線斜率快速的增加，反應出對比劑從血循環系統傳送到腎臟 的情況；第二階段，變化減緩且呈線性(linear)增加，最後會達到整個 MRU 的 最大值，此過程大概經過二至四分鐘，在這個階段代表較多的對比劑會經由血 液進入到腎臟，而較少的對比劑從腎臟離開到收集系統；第三階段，可以用來 計算單腎臟功能(single-kidney function)，利用基準點判斷下降的速度可反應出 對比劑從腎臟到離開收集系統的變化速度。 基於右腎比左腎的比例來表示腎功能，會考慮的兩個參數：(a)第二階段的 斜率(b)組織的功能體積；會考慮第二個參數的原因是：MRU 的腎圖主要是用 腎臟上灰階值的平均量來表示的。 在 MRU 中，腎功能的計算根據在右腎跟左腎第二階段斜率的比值，此部 份只考慮腎每單位所提供的功能而不是全部的加總。於是在 MRU 的腎圖中，
我們將計算第二階段中曲線下面的面積。通常在右左兩邊的腎圖中，第二階段 的起點與終點應該是同步，但是可能有誤差的問題，如當兩邊終點不同，也就 是達到灰階值最高點時間不同，我們則取先達到最高點灰階值的那個時間作為 終點；分別計算出左右腎的第二階段下面的面積再分別乘上左右兩邊腎臟組織 的面積，每邊的比例為自己算出來的積除以兩邊積的加總，此即為單腎臟功能。 圖 2-3 腎功能計算示意圖
2.2
傅立葉轉換(Fourier Transform, FT)
2.2.1
離散傅立葉轉換(Discrete Fourier Transform, DFT)
傅立葉分析是信號分析及處理的一個工具，因為在影像上信號是離散的形 式而非連續，故使用離散傅立葉去處理，以下為離散傅立葉轉換的式子： 對於 N 點序列{f[x]}0≤x<M ，它的傅立葉轉換為
∑
− = − = 1 0 2 ] [ 1 ] [ M n ux M i x f M u Fe
π , u = 0, 1,…, M-1. (1)其中 e 為自然對數的底數，i 為虛數單位。 其反離散傅立葉轉換為：
∑
− = ∧ = 1 2 [ ] ] [ N nk N i n k xe
x
π 0 n , u = 0, 1,…, M-1. (2)2.2.2
快速傅立葉轉換(Fast Fourier Transform, FFT)
在實際應用上，我們選擇由Cooley及Tukey[1]所提出的快速傅立葉轉換
(Fast Fourier Transform, FFT)演算法來做運算，將原來時間複雜度為O(N2)提升
至O(NlogN)，是基於一種叫做連續加倍法(successive doubling method)所推演出 來，為了簡化(1)表示式，則：
∑
− = = 1 0 ] [ 1 ] [ M n ux M f x M u FW
(3) 其中W
M =e
−i2π/M 且假設M =2
n，其中 n 為正整數，M 又可表示為： M = 2K 其中 K 亦為正整數，故取代式子(3)，為：∑
− = = 2 1 0 2 [ ] 2 1 ] [ K n ux K f x K u FW
= [1 [2 ] 1 [2 1]] 2 1 1 0 ) 1 2 ( 2 1 0 ) 2 ( 2∑
∑
− = + − = + + K n x u K K n x u K f x K x f KW
W
(4) 然而，從(3)得知W
22uxK =W
uxK ，(4)可重寫為： = ] [u F [1 [2 ] 1 [2 1]] 2 1 1 0 2 1 0 2∑
∑
− = − = + + K n u K ux K K n ux K f x K x f KW
W
W
(5) 定義： Feven[u]=∑
− = 1 0 ) 2 ( 2 [2 ] 1 K n x u K f x KW
, u = 0, 1,…, K - 1 及 (6) Fodd[u] =∑
− = + + 1 0 ) 1 2 ( 2 [2 1] 1 K n x u K f x KW
, u = 0, 1,…, K – 1 (7) 簡化(5)式子為：2[ ( ) ( ) ] 1 ] [ 2
W
F
F
even u odd u uK u F = + (8) 又因W
uM+M =W
uM及W
u2M+M =−W
u2M，則 ] ) ( ) ( [ 2 1 ] [ 2W
F
F
even u odd u uK K u F + = − (9)2.3
高斯低通濾波器(Gaussian Lowpass Filters, GLPFs)
當影像轉換到頻率領域時，我們可以針對頻率的不同作相對的處理；通常若 是物體本身或背景則會處於低頻部份，若是雜訊或邊界則屬於高頻；而高斯低通 濾波器是讓低頻的部份通過並消除高頻，且比起理想低通濾波器(Ideal Lowpass Filters, ILPF)更為平滑自然；因我們只使用到一維，故僅列出它的一維的轉換函 數為
e
D u D u H 0 ) ( ) ( 2 2 2 − = (10) D(u)為離中心頻率的距離，D0為可調整之參數。 圖 2-4 高斯低通濾波器圖2.4
Canny 邊緣偵測(Edge Detection)
找出圖像的邊緣一直以來是基礎且重要的影像處理課程之一，而邊緣在影像 的某區域中是處於高對比的部份，我們可以利用此特性定義出邊緣偵測器，將一 幅影像過多及比較不重要的資訊過濾，留下影像中未來可以幫助分析的邊緣。
於 1986 年 Canny[2]提出邊緣偵測濾波器的評估準則，它是首先被提出明顯 低錯誤率的邊緣偵測濾波器，以下為提出的三種評估準則： 1. 正確的邊緣偵測能力 - 避免沒找到真正的邊緣以及避免標示雜訊為邊緣。 2. 好的定位(localization)方式 – 被判定的邊緣點應儘可能為真實邊緣的中心位 置。 3. 對於單一邊緣只有唯一的響度(response) – 儘可能減少邊緣附近區域最大 (local maxima)的數量。 根據以上三種準則發展出一套 Canny 邊緣演算法步驟： I. 一開始我們會利用濾波器過濾掉影像上的雜訊，通常使用高斯濾波器 scale σ對影像執行迴旋積(convolution)，假設 f 為影像之函數，G 為二維的高斯 函數，對於影像 f 上的點(x,y)：
e
y x G σ2 2 2 2 + − = (11) II. 而邊緣(x,y)的法向量 n 為 | ) * ( | ) * ( f G f G n ∇ ∇ = (12) III. 使Gn為在n方向G的一階導數 G n n G Gn = ⋅∇ ∂ ∂ = (13)IV. 邊緣點定義為在 n 方向上的 local maximum，則 0 * = ∂ ∂ f Gn n (14) 用(13)式子代入 Gn，則 2 * 0 2 = ∂
∂
n
G f (15) V. 消除雜訊及平滑原始影像之後，接著計算邊緣的強度(strength)，為 | ) * ( | | * |Gn f = ∇ G f (16)VI. 由(12)得知邊緣的法向量方向，垂直法向量方向即可得知邊緣的方向，並可 沿著此方向去尋找下一點的邊緣。
VII. 應用 non-maximal suppression 方法找出邊緣的位置，此方法滿足以上評估準 則 1.及 2.，為： a. 對每一點 C(x, y)，選定垂直於本身方向兩個側邊的鄰近點，為 A、 B； b. 若 M(A) > M(C)或 M(B) > M(C)，則 C 不為邊緣，且設 M(C) = 0， M 為強度函數； c. 輸出強度不為 0 的邊緣。
VIII. 應用 hysteresis thresholding 方法消除毛邊並接續邊緣的斷點 a. 定義兩個thresholds，Thigh和Tlow；
1. 若像素(x,y) > Thigh，該像素標示為strong； 2. 若像素(x,y) <= Tlow，該像素標示為weak；
3. 剩下其他像素接標示為 candidate。
b. 若像素被標為 weak，則略去；若被標為 strong，則輸出為邊緣像素。 c. 若像素被標為 candidate，則判斷判斷是否與 local maxima 相連的邊
緣方向穿過，若是，則輸出為邊緣像素。 d. 若像素被標為 candidate，判斷是否與 strong 像素相連接，若是，則 輸出為邊緣像素。
2.5
影像匹配(Image Matching)
2.5.1 主要演算法為： 1. 建構遮罩(mask)：在資料中選取腎臟灰階值對比最高的影像，畫出腎臟的 邊緣，利用這張影像產生二元(binary)圖，在腎臟邊緣內的像素給予 1 值， 而在外邊則給 0；這張二元圖則為這組資料的遮罩影像；這張遮罩用來跟 不同時間的影像作像素對應相乘，相乘之後的結果叫做 masked images。2. 前處理：在選取的部份事前用 filter 處理，強化影像的資訊，增加最後結 果的正確性。
3. 計算相似值(similarity)：演算法有幾個步驟，每次都會用參考影像(reference
image)與搜尋影像(search image)作比較，計算搜尋影像(search image)在不 同移動之後所得的相似值，最高的相似值所對應的移動為我們估計搜尋影 像對參考影像所移動的值，且記錄此估計值；接著將原本搜尋影像改成參 考影像，並且根據剛剛估計的移動值將遮罩移動到目前的參考影像上，再 對它的下一張影像計算相似值，依此類推，計算到最後一張影像為止。 4. 調整影像(shift)：將求得每個移動值對相對應的影像作校正。 2.5.2 計算相似值的公式推導： 1. Template matching：在一張影像中有一些不同的物體，現在我們有上面物體 的其中一個物體，稱作 template，我們想搜尋它在影像上的位置，通常我們 利用影像與 template 移動之後的對應點灰階值差去計算誤差測量 D(m,n)，且 -M m≤ ≤M，-N n N，m、n 為 template 對影像的偏移量；通常誤差測量是 用平方差(mean-square difference)所定義，可表示為： ≤ ≤
∑∑
−
−
−
= j kn
k
m
j
T
k
j
F
n m D( , )[
(
,
)
(
,
)]
2 (17) 在(17)中，F(j,k)為被搜尋的影像，T(j,k)為 template；且又可表示成； ) , ( ) , ( 2 ) , ( ) , (m n 1 m n 2 m n 3 m n D =D
−D
+D
(18) 在(18)中，∑∑
= j kk
j
F
D
1(m,n)[
(
,
)]
2 (19)∑∑
− − = j k n k m j T k j F n mD
2( , ) ( , ) ( , ) (20)∑∑
−
−
= j kn
k
m
j
T
D
3(m,n)[
(
,
)]
2 (21) ) , ( 3 m nD
代表 template 能量的總和而且它的值不會因為(m,n)的不同而改變， 而 會因 template 位置不同導致求出的能量總和不同，在(4)這項我們可視為影像與 template 之間的交互關係(cross correlation)。
) , ( 1 m n
D
因為我們希望誤差測量越小，故期望D2(m,n)值要大而D1(m,n)值小，根據上面的推斷，我們可以把式子改寫成求他們之間的交互關係，為
∑∑
∑∑
− − = = j k j k FTk
j
F
D
D
R
n k m j T k j F n m n m n m)]
,
(
[
2 1 2 ) , ( ) , ( ) , ( ) , ( ) , ( (22) 改寫之後，我們期望算出來的值是最大。 2. Image registration：由於我們現在所做的並不是 template matching，而是 image
registration，template 變成了另外一張跟原本大小一樣的影像，在做偏移時， template 的能量並不是相同的，所以此時必須將它考慮進去；此外，影像有 時會因拍攝光線過強造成影像震幅(image amplitude)變大，而導致錯誤的估計 值，故我們將欲比對的影像及比對的影像的向量單位化，產生像餘弦的交互 係數(cosine-like correlation coefficient)，亦稱相似值(similarity)：
]
)
,
(
[
]
)
,
(
[
1 2 2 2 ) , ( ) , ( ) , ( 1/2 1/2 2 1 ∑ ∑ ∑ ∑ − − =−
−
∑∑
j k j k n k m j k j n m Rn
n
m
j
F
k
j
F
F
F
j k (23) 在我們應用上，會選取影像上想比對的部份，所以要考慮遮罩的範圍，於是在 式子(23)中加上 M(j,k)，改寫為：]
)
,
(
)
,
(
[
]
)
,
(
)
,
(
[
1 2 2 2 ) , ( ) , ( ) , ( ) , ( 1/2 1/2 2 1 ∑ ∑ ∑ ∑ − − =−
−
∑∑
j k j k n k m j k j k j M n m Rn
k
m
j
F
k
j
M
k
j
F
k
j
M
F
F
j k (24)2.6
次影像(Sub-pixel)
本實驗資料間的位移通常很小，所計算出來的值常在一至兩個像素之間， 但由結果比較原始資料可看出有時有校正過度的情況，於是我們將原本資料的解 析度 256×256 利用雙線性內插(bilinear interpolation)的方式提昇至 512×512；以下 為各位置內插的公式： 1. 欲計算的像素在原始影像的內部：圖 2-5 次像素計算示意圖 1 E = 1 / 4 (A + Ｂ + C + D) (25) 2. 欲計算的像素在原始影像的左邊界或右邊界 圖 2-6 次像素計算示意圖 2 C = 1 / 2 (A + B) (26) 3. 欲計算的像素在原始影像的上邊界或下邊界 圖 2-7 次像素計算示意圖 3 C = 1 / 2 (A + Ｂ) (27)
三、實驗流程
在本章中我們會利用第二章已經提及到的方法，建構兩種不同的演算法，第一種方 法我們將會在頻率領域上作處理，第二種方法則是應用在空間領域上，實驗的資料是腎 功能有問題的 MRU，根據前一章計算腎功能的公式，最後比較這兩種演算法算出來的 差異，並且比較它們之間優缺點。 此處實驗資料有五組，以下為資料內容的說明： 資料編號 維度 說明 (1) 256×256×170 單邊水腎 (2) 256×256×170 單邊水腎 (3) 256×256×200 雙邊水腎 (4) 256×256×160 雙邊水腎 (5) 256×256×170 單邊水腎 表 3-1 實驗資料(1) (2)
(3) (4)
(5)
3.1
第一種演算法
3.1.1流程圖
首先我們先畫出這個演算法的實驗的步驟，再針對實驗過程中所產生的結 果及問題做詳細的說明： 開始 快速傅立葉 轉換 高斯濾波器 反快速傅立 葉轉換 結束 圖 3-2 演算法一流程圖 3.1.2 針對時間軸做快速傅立葉轉換原始資料：
圖 3-3 原始資料示意圖
經過 3D 對 z 軸的快速傅立葉轉換：
圖 3-4 經高斯低通濾波器
3.1.3 利用高斯濾波器讓資料變化平滑
(1) 參數為 10
(2) 參數為 20
(3) 參數為 30
(4) 參數為 40 圖 3-8 高斯參數 40 反傅立葉之結果 當高斯參數越小時，沿著時間軸的影像表現的越平滑，原本因呼吸而上下移動 的現象會被改善，但從另一個角度來看，越多原始資料的資訊會被過濾，可能會增 加誤差值，所以我們希望既能改善移動亦能保留大部分的資訊，因此我們測試了四 個高斯參數。 以下為五組資料，各分為左右兩邊腎臟，為皮質部位在不同時間的灰階值變 化，橫軸座標代表影像拍攝次序的編號，每一張影像的間隔時間為 10 秒鐘，縱軸
為皮質部位所有像素的灰階值平均： (1) 左 右 (2) 左 右 (3) 左 右 (4) 左 右
(5) 左 右 圖 3-9 MRU 曲線之結果 根據以上資料的顯示，他們的高斯參數皆會在 40 的時候變化驟減趨近於 收斂，因此之後我們在處理以後的資料時就使用這個參數值。 另外，在每組 MRU 中，可看出一開始的平均灰階值會偏高而最後的部份 則會偏低，主要的原因是運算的時候是對資料以環狀方式做處理，導致前後有 明顯不平滑的情狀。
3.2
第二種演算法
3.2.1流程圖
首先我們先畫出這個演算法的實驗的步驟，接著再針對實驗過程中所產生 的結果及問題作詳細說明：建構遮罩 資料前處理 初始參考資料 計算兩影像間的相似值 紀錄最高相似值的位移 指定參考影像前一張 影像為搜尋影像 指定參考影像後一張 影像為搜尋影像 計算兩影像間的相似值 紀錄最高相似值的位移 搜尋影像為第 一張影像 搜尋影像為最 後影像 搜尋影像改成 參考影像 搜尋影像改 成參考影像 調整全部 影像 開始 結束 圖 3-10 演算法二流程圖 3.2.2 建構遮罩
遮罩的功用主要是縮小考慮範圍，增加它的正確性，我們利用人工的方式圈 選需要校正的部份，且必須包含到全部影像的右(左)腎，所以在左右腎我們都會 各建構一個遮罩： 圖 3-11 左右腎之遮罩圖 3.2.3 資料前處理 (1) canny edge 由於計算相似值主要是根據兩張影像間的灰階值來比對，但拍攝相鄰兩 張影像並非相同的 sensor，相同處的灰階值幾乎都不相同，會造成比對錯誤 的情況，為了避免這樣的結果，我們將原始影像留下重要的資訊並簡化他們 的灰階值，於是利用 canny edge detector 將邊界找出來並讓它的灰階值為 1 其餘則為 0。
圖 3-12 前處理後
(2) 次像素
經過第一次的前處理之後，校正的結果很合理也有改善，但有些地方會 有校正過度的情況，如：
圖 3-13 有無前處理及次像素之比較圖 這是其中連續兩張影像的合成，前一張為紅色，後一張為綠色；觀察藍 色視窗的部份，從(a)來看，後一張影像往上偏移，經過 canny 方法的校正之 後，後一張影像有偏下的情況，若改加上次像素的方法，有修正半個像素卻 又無矯正過度的情況。 從上來看，為了讓結果更好更精確，我們加上使用次像素的方式，先讓 原始影像放大再去計算相似值，計算之後再將大小縮回原尺寸，即使有移動 半個像素的情況也能計算出來。
3.3
兩種演算法之比較
以上介紹了兩種演算法，顯然的第一種方法過程比較簡單，但因為有使用 到高斯濾波器，它計算出來的數據必定有些誤差，於是我們希望從第二種演算 法來驗證究竟可不可行。 3.3.1 MR urographic 之比較 首先，我們先從五組資料的 MR urographic 來比較，主要衡量的標準就是以 第二階段的面積所求出來的左右腎的比例，因為未來需要使用在醫學統計上面。(1) 20979490 甲、右腎的比較 圖 3-14 資料一之右腎比較 乙、左腎的比較 圖 3-15 資料一之左腎比較 丙、兩腎比例的比較 腎功能比例(%) 演算法 右 左 1 11.35027 88.64973 2 11.46519 88.53481 絕對值差 0.11492 0.11492 表 3-2 資料一兩腎比例之比較
(2) 20505907-20060404 甲、右腎的比較 圖 3-16 資料二之右腎比較 乙、左腎的比較 圖 3-17 資料二之左腎比較 丙、兩腎比例的比較 腎功能比例(%) 演算法 右 左 1 40.21057 59.78943 2 38.62115 61.37885 絕對值差 1.58942 1.58942 表 3-3 資料二兩腎比例之比較
(3) 20588343-20060404 甲、右腎的比較 圖 3-18 資料三之右腎比較 乙、左腎的比較 圖 3-19 資料三之左腎比較 丙、兩腎比例的比較 腎功能比例(%) 演算法 右 左 1 60.34557 39.65443 2 59.94071 40.05929 絕對值差 0.40486 0.40486 表 3-4 資料三兩腎比例之比較
(4) 20905437-20060412 甲、右腎的比較 圖 3-20 資料四之右腎比較 乙、左腎的比較 圖 3-21 資料四之左腎比較 丙、兩腎比例的比較 腎功能比例(%) 演算法 右 左 1 59.79936 40.20064 2 58.7364 41.2636 絕對值差 1.06296 1.06296 表 3-5 資料四兩腎比例之比較
(5) 20925363-20061108 甲、右腎的比較 圖 3-22 資料四之左腎比較 乙、左腎的比較 圖 3-23 資料五之左腎比較 丙、兩腎比例的比較 腎功能比例(%) 演算法 右 左 1 67.77422 32.22578 2 68.74374 31.25626 絕對值差 0.96952 0.96952 表 3-6 資料五兩腎比例之比較
在分別比較左右腎比例時，兩演算法所呈現出來的結果看起來都很相近，演算 法一因為用了濾波器過濾掉不平滑的資訊，所以曲線較平滑；而演算法二用的是原 始資訊，且相鄰影像的灰階值不同，致使曲線皆有上下起伏的情況。 從上面結果來觀察，它們雖然很相近但終究有誤差，更進一步的，我們將它們 的結果數據化，計算出兩腎之間的百分比，它們的絕對值差皆在 2%以內，從這裡 我們約略可說它們的結果幾乎相近。 3.3.2 時間之比較 (1) 演算法一時間 = FFT 時間 + 高斯濾波器時間 + IFFT 時間
(2) 演算法二時間 = 選擇 ROI 時間 + 事前處理(Canny edge + Subpixel) 時間 + 影像匹配時間 時間(秒) 資料 演算法一 演算法二 (1) 25.23 1644.16 (2) 23.3 1615.24 (3) 25.5 1990.55 (4) 20 1588 (5) 24.49 1708.75 平均 23.704 1709.34 表 3-7 兩演算法時間之比較 演算法二的過程較繁雜而且它所花的時間也比較長，幾乎是演算法一的七十倍 以上，所以我們之後要處理相同問題時，將利用演算法一來處理資料。
四、實作與結果
本資料來自於長庚兒童醫院，實驗共六十二組資料，有五十四組資料屬於單 邊水腎，剩餘八組資料為雙邊水腎。
4.1
實驗環境
(1) Multi CPU：Intel(R) Core(TM)2 CPU 6300(1.86GHz)
Intel(R) Core(TM)2 CPU 6300(1.86GHz)
(2) RAM：2 GB
(3) 作業系統：Microsoft Windows XP Professional
4.2
MRU 電腦協助計算與放射科醫師結果比較
4.2.1 放射科醫師之作法 與電腦協助計算的最主要的分別是： (1) 實驗的資料沒有事前校正，會隨著呼吸而上下移動。 (2) 由放射科醫師圈選皮質，如圖 4-1。 (3) 在計算各別腎功能時，沒有精確的計算而是以肉眼判斷最大斜率。 (4) 計算 MR urographic 的第二階段下的面積是以計算三角型面積去估計，如圖 4-2。 (5) 使用的是訊號值而非灰階值。圖 4-1 放射科醫師圈選之皮質
圖 4-2 放射科醫師使用之 MRU 4.2.2 比較
圖 4-3 電腦與放射科醫師結果之比較 誤差比例% 0 至未滿 5 5 以上至未滿 10 10 以上至未滿 15 15 以上 組 49 9 4 0 表 4-1 誤差比例之分配表 4.2.3 討論 上面有提到五個放射科醫師與電腦協助計算不同之處，由(1)來看，雖然放 射科醫師沒有校正，但大多數資料並不會上下移動的很劇烈，而比例相差超過 11 的資料也沒有這種情形；而在(2)的部份，電腦協助計算的部份為先取 threshold 再刪去不要的部份，此電腦協助操作比較節省時間但因人為的操作還是存在誤 差，在這些例子中，是影響結果最大的因素；在(3)中，很有可能以肉眼判斷因 相鄰的點而產生 5 以上的差距；在(4)中，用三角型面積公式去估計曲線下面的 面積也有可能產生誤差；除此之外，在(5)中放射科醫師使用的訊號值，電腦協 助計算使用著則是原始資料 2byte，也會有部份的誤差；因在利用 MRI 去計算腎 功能還沒有辦法完全去除人為因素，所以我們希望盡可能減少人為的介入以接近
真實的腎功能。
4.3
MRU 與核子醫學結果之比較
目前醫學上大多數醫師還是採取核子醫學的結果作為主要參考值，但是在我 們實驗中，很多資料會出現顯而易見的錯誤，如圖 4-4，左腎皮質有功能的地方 明顯高於左腎，但在核子醫學方式求得的卻是：右腎比左腎的比值為 52：48， 在 MRU 的部份則為 37.55:62.45，雖然我們無法確定我們求出的值是否為真實的 值，但從圖中觀察卻比較合理，於是我們希望能從比較 MRU 與核子醫學的結果， 整理出核子醫學可能出錯的情況。 圖 4-4 核醫錯誤之範例 (1) 兩結果的相關圖圖 4-5 MRU 與核醫結果之比較 上頁圖 4-5 列出有水腎或水腎較嚴重的那邊算出的值，根據圖 4-5 的 統計結果，高達約 47 組資料核醫的值較高，整體來看核醫算出來的值估 得比較高，其中數列 2 共有 16 組數則是錯誤的值，因為目前無法得知標 準值，只能從明顯的錯誤去判斷結果是不是錯誤的，所以我們只能說 62 組中至少有 16 組是錯的，且集中於 40 至 50 之間。所以當兩邊腎功能相 近時，核醫容易產生錯誤。 (2) 兩腎長度差與兩方法之差 由於正常腎臟的長度可以代表年齡，故我們假設水腎長度或較嚴重水 腎的長度為 R，另一則為 L，兩腎臟長度差為(R-L)/L ，如圖 4-6， 數列 1 為單邊水腎，且核醫的值小於 50%，數列 2 為核醫所求的值超過 50%，數列 3 為雙邊水腎，且核醫的值小於 50%，而數列 4 為核醫所求的 超過 50%。 % 100 ×
從圖 4-6 可得知，長度差在 40%為一個臨界點，以此臨界點我們將 MRU 與核醫所求出來的值分為兩個 Group，如圖 4-7 中得知，長度差在 40%以內時，核醫的結果幾乎大於 MRU，且數據分布較集中；如圖 4-8 可知，長度差超過 40%之後，也就是水腎腫脹厲害時，核醫的方法會估得 較高或較低，數據相當不集中，此時 MRU 之方法須從中取出過薄的皮質， 也不好判斷，其實兩種作法都不是相當可靠，所以當水腎嚴重時，可能需 要參考兩種方法所求出來的值。 圖 4-6 MRU 與核醫相差值與腎長度之相關圖
圖 4-7 長度差(%)低於 40 之圖
圖 4-8 長度差(%)高於 40 之圖
另外，在圖 4-9 中，左腎其實還有微薄的皮質，並有功能存在， 但核醫無法偵測到其腎功能，但 MRU 還能偵測到其腎功能，在某些時候，
MRU 求出的值更具可信度。 圖 4-9 核醫偵測不到腎功能之範例 (3) 水腎中皮質佔全腎臟之比例與兩方法之差 上段我們以長度的方法代表水腎的嚴重程度，因為目前並無水腎嚴重程 度統一的說法，故我們接著以水腎中皮質面積佔全部腎臟面積的比例再次去 做統計，若比例越小代表水腎嚴重度大，也就是成反比。 如圖 4-10 所示，在 50%/的時候為一個臨界點，當在 50%以上的時候， 核醫所求的值幾乎是大於 MRU，而當百分比愈小的時候，這兩種方法的差異 性愈大；不過這個統計圖並不能有很大的意義，因為它的數據過於分散，隱 約呈現於一個金字塔，故我們只能這樣推測它的結果。
