抗核抗體免疫螢光顯影分析之處理方法評估

(1)

抗核抗體免疫螢光顯影分析之處理方法評估

葉進儀* 李貞儀

國立嘉義大學資訊管理學系

摘要

免疫風濕科在檢查自身免疫疾病的主要依據是抗核抗體 (antinuclear autoantibodies, ANA)，醫師會藉由人類表皮細胞癌 (human epithelioma type 2,

HEp-2) 之細胞培養，透過免疫螢光顯影技術來辨識 ANA，以預測病人的疾

病。目前辨別樣式是仰賴專家觀察螢光顯微鏡下的載玻片來完成，此方式需要具有高度專業之技術人員來操作且耗時，因此自動化分析免疫螢光顯影影像有其必要性。目前研究自動化分類 ANA 螢光影像流程的技術眾多，本研究從眾多方法中找出分類結果正確率較高之組合方法，進而協助醫師進行疾病之診斷。實驗結果顯示使用 Canny 邊緣檢測，經紋理分析得到特徵，由支援向量機遞迴特徵消去法進行特徵選擇，最後結合 SVM 分類方法具有最高的正確率 (97.00%)。

關鍵詞：HEp-2 細胞，邊緣檢測，紋理分析，特徵選擇，分類。

EVALUATION OF PROCESSING METHODS FOR ANALYZING ANA IMMUNOFLUORESCENCE IMAGES

Jinn-Yi Yeh* Chen-Yi Lee

Department of Management Information Systems National Chiayi University

Chiayi, Taiwan 600, R.O.C.

Key Words: HEp-2 cells, edge detection, texture analysis, feature selection, classification.

ABSTRACT

Experts of immunology and rheumatology department inspect auto- immune disease by recognizing antinuclear autoantibody (ANA) patterns.

Physicians diagnose patients’ disease by inspecting HEp-2 cells for identification of ANA. Up to now, identification of ANA is completed by inspecting slides using a fluorescent microscope. This manual procedure requires highly specialized technicians and is very time-consuming.

Therefore, it is necessary to have methods for classifying ANA images automatically. This research attempts to find the cascaded method with the highest accuracy from among different processing methods. These methods include edge detection, feature selection, and classification. The experimental results show that the optimal cascade methods are the com-

*通訊作者：葉進儀，e-mail: jyeh@mail.ncyu.edu.tw

Corresponding author: Jinn-Yi Yeh, e-mail: jyeh@mail.ncyu.edu.tw

(2)

46 技術學刊第二十八卷第一期民國一○二年

bination of Canny and Support Vector Machines (SVM). The accuracy rate is about 97.00%.

一、前言

免疫風濕科主要診斷因身體免疫系統失調或功能異常所造成之疾病症狀，這些疾病可能影響身體所有的器官，表現出來的症狀也相當多樣化，常見的疾病有紅斑性狼瘡、硬皮症、類風濕關節炎、僵直性脊椎炎等，而抗核抗體 (antinuclear antibody, ANA) 是血中的抗體對抗細胞核內的抗原，意即對抗自己細胞核內抗原的自身抗體

(autoantibody)。這些自身抗體的存在，與各種不同的免疫

疾病息息相關，對於診斷不同的風濕性疾病具有非常重要的意義。ANA 判讀包括使用螢光顯影 (fluorescence) 技術之螢光強度指數 (intensity factor, IF) 亦稱效價 (titre) 及樣式 (pattern) 兩大要素。Boland 等[1]首先提出將螢光細胞影像利用 Zernike 多項式計算與 Harlick 材質分析來產生圖形樣本的特徵值，再以統計和倒傳遞類神經網路 (BPNN) 比較五種次細胞或蛋白質結構分類。Murphy 等[2]將過去的研究加上類神經網路來學習分辨，進一步探討類神經網路隱藏層的節點與階層數對於螢光細胞影像辨識準確率的比較，另外除了計算先前兩種特徵值：Zernike 矩特徵 (Zernike moments)、Haralick 紋理特徵 (Haralick texture

features)，也加入了邊緣檢測與外型輪廓幾何的特徵計算，

增加樣本影像相似時的可辨性。

利用此免疫螢光顯影的技術於人類表皮細胞癌 (human epithelioma type 2, HEp-2) 上，是目前免疫風濕科在辨別 ANA 類別最主要的方式。最早使用免疫螢光顯影技術於辨識 ANA 是在 1957 年，此技術是利用 HEp-2 細胞來表示主要的基質，具有以下特性：(1) 非常均勻的細胞樣式；(2) 細胞本身具有很大的細胞核；(3) 細胞具有代表性的細胞質；(4) 在任何細胞週期的狀況下都會產生此細胞[3]。因為具有這些特性，所以容易辨別出 ANA 的類別，

如果給予一百種不同的 ANA，觀察 HEp-2，即可辨識出超過三十種不同的細胞核及細胞質樣式。然而目前處理的流程，需仰賴高度專業的技師或醫師，以人工的方式來觀察螢光顯微鏡下的載玻片後，以主觀經驗來判斷 ANA 的類別，此處理流程有幾個缺點：(1) 因需主觀判斷類別，所以需要高度專業的人員才能進行；(2) 缺乏全自動化的解決辦法，以人工辨識相當耗時；(3) 缺乏可靠的標準化文件[4]。

由於有上述缺點，所以近年來有諸位學者致力於將此流程自動化，例如 Hahm and Anderer 等[5]探討了在對 ANA 螢光顯影做影像處理前的準備程序，實驗了使用不同的凝固、染色藥劑對保存細胞樣式的差異，期望找出哪個凝固程序在保存跟疾病有關的抗原上有最佳的結果，因而最適用在自身抗體的抗核抗體測試上。最早對 ANA 螢光顯影

表一分類方法優缺點

分類方法優點缺點

類神經網路 1. 平行處理 2. 錯誤的容忍度高 3. 聯想記憶高 4. 處理連續型資料較

強

5. 解決最佳化問題 6. 應用性較廣

1. 產生之法則較不易解釋維護 2. 有無限多組解 3. 神經元個數太多會

影響計算正確性 4. 耗費大量記憶體決策樹 1. 分類出的規則容易

理解

2. 資料運算處理時間不長

3. 可處理類別型字串資料

1. 處理連續型字串資料較弱

2. 處理時間序列資料需做離散化 3. 資料類別太多時錯樸素貝氏結合事前機率與樣本誤率高

機率，可得到更多資訊假設所有屬性彼此獨立，實際不可能支援向量機 1. 學習速度快

2. 具良好統計學基礎 3. 處理連續型資料較

強

4. 可處理高維度資料

1. 不易建立最佳 Kernel function 2. 不易解釋分類過程

做分類的是 Perner [6]，使用之特徵為統計型特徵 (灰階平均、細胞周長、平均面積等)，使用之分類器為決策樹分類器，分類錯誤率為 13.33%。黃嘉政[7]使用 Canny 邊緣檢測得到細胞輪廓，得到輪廓後使用紋理分析萃取出影像特徵，其特徵亦為統計型，再輸入至類神經網路做分類；然而其影像缺乏對比，影響了邊緣檢測，分類錯誤率為 15.62%。Huang 等[8]使用分水嶺切割 HEp-2 細胞，分類結果之敏感度 (有疾病者被預測為有疾病者之機率) 為 94.7%。

Sack 等[4]實作一套裝軟體來辨識 HEp-2 細胞樣式，

套裝軟體的基礎是新的影像分析及特徵選擇演算法，影像從相機輸入轉換成灰階影像，細胞從影像中被切出，描述這些不同的單一細胞的特徵被存在資料庫，對產生的資料集做機器學習找出相關的特徵，以便能夠正確的分類自身抗體，最後用學習演算法進行特徵選擇來做分類器，以軟體分類之錯誤率為 16.91%。

就影像處理而言，學者大多使用的邊緣檢測方法有

Sobel、高斯的拉普拉斯、LoG、Canny、小波轉換等；形

態分析有膨脹、侵蝕、斷開、閉合等；紋理分析有結構型、

模式型、統計型、轉換型等四種不同之特徵參數[9]。就分類流程而言，分類前所使用的特徵選擇有多種機制，例如：

Information Gain、ReliefF；各學者選用的分類器也不同，

例如：Nattkemper 等[10]以類神經網路將 HEp-2 細胞分類，

Perner等[11]則是以決策樹將 HEp-2 細胞分類，各分類方

法皆有其優缺點 (表一)。

(3)

圖1 研究流程圖

由此可見，自動化流程的各階段均有多種方法，且各研究所使用之影像亦不同，所以無從比較何種方法具較佳分類結果。因此，本研究評估邊緣檢測器、紋理分析、特徵選擇方法及分類器，藉以找出分類結果最正確之組合方法，輔助醫師達到加速診斷疾病的功能。

二、ANA螢光影像分析方法

本研究擷取化學滴定劑檢查所得之病人血液影像，進行處理、分析，其架構分為五個階段 (見圖 1)：(1) 影像前處理，將影像中之灰階對比增強；(2) 針對 ANA 影像，使用 Canny、小波轉換邊緣檢測來圈選細胞質與細胞核之輪廓，擷取細胞所在位置。擷取出細胞位置後再透過形態分析，得到更加完整的邊緣；(3) 將所圈選出之細胞進行紋理分析，擷取細胞紋理特徵；(4) 經過紋理分析之後，從特徵集合中選出重要特徵，將細胞進行分類；(5) 進行分類效能評估，找出最佳之抗核抗體影像分類系統。以下針對影像前處理、邊緣檢測、紋理分析、特徵選擇、分類、

評估逐一描述：

1. 影像前處理

在許多情況下都可能產生對比度低的影像，例如照明不夠、影像感應器動態範圍不夠，或者在影像攝取期間透鏡光圍的錯誤設定等。「對比度擴展 (contrast stretching)」

是增加待處理影像的灰階範圍，將高於 m 的灰度變亮、低 於 m 的灰度變暗 (圖 2)，圖 2(b)所示為極限情況，將高於 m的灰度變白，低於 m 的灰度變黑，轉換成二值影像。

2. 邊緣檢測

雜訊的去除與細胞邊緣的保留是互相衝突的，如果使

圖2 用於對比度擴展的灰階轉換函數[12]

用低通濾波器去除雜訊，會使得影像變得非常模糊，影響到邊緣檢測的位置與準確度；而利用高斯濾波器減少雜訊會產生邊緣誤差，必須在去除雜訊與邊緣檢測衝突中找到平衡。目前廣為人知的邊緣檢測運算子有 Sobel、LoG、

Canny [9]，其中以 Canny 最常被使用且檢測結果最好，因

為 Canny 藉由尋求像素之梯度的最大值來找邊緣，梯度的計算是用高斯濾波器的導數，此方法使用兩個臨界值來檢測強邊緣與弱邊緣，並只輸出與強邊緣有相連接的弱邊緣，因此此方法較有可能檢測出真正的弱邊緣。Li and

Zou [13]提出以小波轉換來做邊緣檢測之結果與 Canny 相

仿，因此本研究所要評估的邊緣檢測方法是 Canny 及一階二維小波轉換。

(一) Canny 邊緣檢測

Canny的邊緣檢測能夠將影像去除雜訊使得平滑化，

使用時須選擇適當的高斯濾波器，其中標準差σ可以控制濾波器的範圍，σ 越大抑制雜訊的能力越強，但 相對地資訊也可能因此而被過濾掉，反之σ 越小雖可獲得較多的邊緣資訊但相對地不必要的雜訊也可能被保留下來，所以需要謹慎選擇σ值[14]，以下說明此演算法：

(1) 高斯濾波器 (gaussian filter) 主要的目的就是將影像雜訊和不必要的細節和紋理去除掉，之後再找出影像中的邊緣。

(2) 將每一個點都計算出局部梯度 g(x, y) 與邊緣方向

θ^{(x, y)}，然後邊緣點被定義為在此梯度方向上有局

部最大強度的點。

(3) 非最大值的抑制 (non-maximal suppression) 為在梯度影像中，對於每一個被找到的邊緣點作為梯度方向的搜尋，與其方向上所有的像素點做比較，如果某像素為該梯度方向上所有梯度值之中的最大值，則將這個像素點標記為邊緣上的候選點，否則標記為非邊緣點；經過此方法處理過後的影像，會得到一個邊緣細線化的結果。

(4) 為了補救邊緣斷裂的問題和刪除不正確的邊緣 點，Canny 給定兩個臨界值 T1和 T2，其中 T2 > T1，以用來增加整個邊緣檢測上的彈性，若邊緣上的候 選點，梯度值大於 T2則判定為邊緣點，對於

(4)

f(n)

2 ↓

2 ↓ h(k)

g(k)

2 ↓

2 ↓ h(k)

g(k) LL

LH

HL

HH 2 ↓

h(k)

g(k)

L

H

圖3 二維離散小波轉換[15]

梯度值介於 T1和 T2之間的點，就會把此點補救回來，

但是若梯度值小於 T1則判定為非邊緣點。

(二) 小波轉換 (wavelet transform) 邊緣檢測

小波轉換對數位信號具有多重解析特性，各頻段間相互獨立，使得局部分析效果良好，其多分辨率能力為影像的邊緣檢測提供了新方法，近年來常被使用在影像壓縮、紋理分析等影像處理上[13]。小波轉換影像通常是以二維平面的方式表現，二維的離散小波轉換流程如圖 3。二維小波轉換是先針對平面影像中水平方向的信號，每一列進行一維的小波轉換，分解成為水平方向上的低頻與高頻信號。然後，再對水平方向高頻和低頻兩部份的垂直部分，每一行依序進行一維的小波轉換。在離散小波轉換時，並不需要實際求出小波函數及比例函數，只要利用高通跟低通濾波器，

經過二分之一的取樣 (down sampling)，就可以求得離散小波轉換。圖 3↓表示取出影像中兩相鄰的像素點的數值，經過運算後得到一個轉換後的數值，h(k) 表示低通濾波器，g(k) 代表高通濾波器。信號經過轉換，分解成為一半表示低頻部份，一半代表高頻的信號。信號還原時，將信號經過高頻或低頻率波的逆轉換，回復成為原始的兩個信號數值。在一階二維的小波轉換裡，影像會被分為 LL、LH、HL、HH，如圖 4。

本研究將離散小波轉換運用在邊緣檢測上，將 ANA 灰階影像經過一階二維離散小波轉換，會得到不同之方向資訊，要結合這些資訊到單一影像的邊緣上，我們將所產生之小波轉換的近似係數取零，計算其反轉換，再取絕對值，即可得到影像邊緣。

3. 紋理分析

本研究從像素資料中取得紋理資訊，由此提供影像的分類，此階段除了應用過去研究所使用的共生矩陣 (co-occurrence matrix)，也考量了遊程矩陣 (run length

matrix) 所計算出之紋理特徵。紋理分析經常應用於影像元

素的灰階值與座標空間中的相對位置，再進行統計方法的演算與分析進一步瞭解影像結構內容的方向性、規則性或粗糙度等。當找出紋理描述後，用來比較不同質地或紋理影像的差異變化。

圖4 一階二維離散小波轉換[15]

共生矩陣的產生方式是統計影像空間內，每一個像素在符合預先設定的條件下（如：水平方向、相隔 1 個像素），出現相同灰階值的次數[3]。本研究共生矩陣採用之角度為：水平方向(0°)、右上對角方向(45°)、垂直方向 (90°)、左上對角方向(135°)以上四個方向；距離為 1、2、

4、6、8 五個向量[16, 17]。例如包含四種不同灰階值 (0~3) 之影像矩陣 I (4 × 4 矩陣)，假設θ = 0°、45°、90°、135°，

且 d = 1 時，將可以萃取出 36 種特徵 (9 種特徵* 4 種不同角度之組合 = 36 種特徵)，因此由共生矩陣產生 180 個紋理特徵。假設影像大小為 M × N，其共生矩陣所包含 的參數 P 是表示一個 n × n 的矩陣，在一張影像中 n 是灰 階值的數量，為了計算上的效率，需要減少灰階值的數量，減少共生矩陣的大小，而矩陣是用來扮演一個累加 器，P[i, j]是計算影像中所有成對的像素，i 和 j 為灰階值 的強度， ( )P i_x =

∑

P i j[ , ]為加總 P[i, j]行的邊際機率矩陣，

( ) [ , ]

P jy =

∑

P i j 為加總 P[i, j]列的邊際機率矩陣。計算公 式如表二所示。

共生矩陣提供了影像中灰階值空間分布情形的資訊，但只有求出共生矩陣是不夠的，必須將它進一步量化為單一數值的形式，方能提供電腦直接分析使用。假如影 像大小為 M * N，本研究將紋理做量化，由下列演算公式 計算分析影像特徵[18-20]。

遊程矩陣是由搜尋整張影像，統計某像素灰階值在給定一方向 (direction) 、距離 (run size) 的條件下計算出來的。本研究採用方向有四：0°、45°、90°、135°，距離為 2 到 7，共四個遊程矩陣，假設 f(i, j)為影像遊程矩陣灰階 i 與遊程距離 j 的數目、G 為灰階強度、R 為最大遊程、nr

為遊程距離總數、np為此影像像素總數、fr(j) 為遊程距離 函數、fg(i) 為影像像素資訊，每一矩陣可計算出十一個紋理特徵，因此由遊程矩陣產生 44 個紋理特徵，如表三。這些特徵為直覺推論而推導出之遊程分佈特性，例如圖 5 顯示正負 0°、45°、90°、135°方向的 8 種權重加總遊程特徵。

其缺點為缺乏理論基礎，且特徵間有高度相關，但應用於分類卻有不錯之績效[21]。

(5)

表二共生矩陣紋理特徵

紋理特徵表示式

平均 (mean, µx, µy)：平均強度的量測。 [ , ]

N

j x

P i j

µ =

∑

N ^M ^{[ , ]}

i y

P i j µ =

∑

M 變異數 (variance, σ²⁾：量測灰階值分散的程度，灰階值分散程度越大，變異數越大。

2 2

( ) [ , ] ( ) [ , ] 2

M N

x y

i j

i u− P i j + −j µ P i j

∑∑

對比度 (contrast)：量測影像中對比的強列程度。 ( )² [ , ]

M N

i j

i−j P i j

∑∑

相關度 (correlation)：量測灰階值線性相依的程度，即是在成對的像素中兩像素之間的關係程度，可顯示紋理的方向性，當多數紋理結構朝某一特定方向排列時，則代表在此一方向的紋理相關度比其它方向高。

( )( ) [ , ]

M N

x y

i j x y

i µ j µ P i j σ σ

− −

∑∑

能量 (energy)：量測紋理的一致性，影像紋理越一致，能量值越大。 [ , ]²

M N

i j

P i j

∑∑

平滑度 (smoothness)：量測一區域中強度的相對平滑度。對於一個常數強度的區

域，R 為 0；而對強度位準上有較大變動的區域，其值逼近 1。 ² 1 1

1 σ

− + 第三矩 (u3)：量度直方圖的歪斜性，對於對稱直方圖，此量測為 0，直方圖向右斜 (對平均) 時為正，而向左斜時則為負。

( , )3 [ , ]

M N

x y

i j

i−µ j−µ P i j

∑∑

均質性 (homogeneity)：當所有灰階值相等 (最均勻) 時，此量測最大並從此開始減少

[ , ]

1 | |

M N

i j

P i j i j

∑∑

+ −

熵 (entropy)：量測隨機性。 [ , ]log [ , ]

M N

i j

P i j P i j

−

∑∑

表三遊程矩陣紋理特徵[21]

紋理特徵表示式

短遊程優勢 (Short Run Emphasis, SRE) ² ²

1 1 1

1 ^G ^R ( , ) 1 ^R _r( )

i j j

r r

f i j f j

SRE=n

∑∑

₌ ₌ j =n

∑

₌ j

長遊程優勢 (Long Run Emphasis, LRE) ² ²

1 1 1

1 1

( , ) * ( ) *

G R R

r

i j j

r r

LRE f i j j f j j

n ₌ ₌ n ₌

=

∑∑

=

∑

灰階非均勻性 (Gray-Level Nonuniformity, GLN) ² ²

1 1 1

1 1

( ( , )) ( )

G R G

g

i j i

r r

GLN f i j f i

n ₌ ₌ n ₌

=

∑ ∑

=

∑

遊程非均勻性 (Run Length Nonuniformity, RLN) ² ²

1 1 1

1 1

( ( , )) ( )

R G R

r

j i j

r r

RLN f i j f i

n ₌ ₌ n ₌

=

∑ ∑

=

∑

遊程百分比(Run Percentage, RP) RP = nr / np

低灰階遊程優勢 (Low Gray-Level Run Emphasis, LGRE) ² ²

1 1 1

1 ^G ^R ( , ) 1 ^G g( )

i j i

r r

f i j f i

LGRE=n

∑∑

₌ ₌ i =n

∑

₌ i 高灰階遊程優勢 (High Gray-Level Run Emphasis, HGRE) ² ²

1 1 1

1 1

( , ) * ( ) *

G R G

g

i j i

r r

HGRE f i j i f i i

n ₌ ₌ n ₌

=

∑∑

=

∑

低灰階短遊程優勢 (Short Run Low Gray-Level Emphasis, SRLGE) ² ²

1 1

1 ( , )

*

G R

i j

r

f i j SRLGE

n ₌ ₌ i j

=

∑∑

高灰階短遊程優勢 (Short Run High Gray-Level Emphasis, SRHGE)

2 2

1 1

1 ^G ^R ( , ) *

i j

r

f i j i SRHGE

n ₌ ₌ j

=

∑∑

低灰階長遊程優勢 Long Run Low Gray-Level Emphasis, LRLGE)

2 2

1 1

1 ^G ^R ( , ) *

i j

r

f i j j LRLGE

n ₌ ₌ i

=

∑∑

高灰階長遊程優勢 (Long Run High Gray-Level Emphasis, LRHGE) ² ²

1 1

1 ( , ) * *

G R

i j

r

HGRE f i j i j

n ₌ ₌

=

∑∑

(6)

S R E

L R E

LGRE SRLGE LRLGE

HGRE SRHGE LRHGE

圖5 八種權重加總遊程特徵[21]

4. 特徵選擇

得到上述特徵後，由於特徵數眾多，會導致分類時學習困難、耗時，且可能存在不相關或相似度高的特徵，所以必須進行特徵數量較佳化之選擇，在特徵評估階段使用以下方法將特徵與予排序：(1) 資訊獲利 (information gain)；(2) 獲利比率 (gain ratio)；(3) ReliefF；(4) 支援向量機遞迴特徵消去法 (recursive feature elimination, RFE)；

(5) OneR。在特徵選擇階段採用統計變異數分析之事後檢

定法加上績效考量，分別選出 10, 20, …, 70 及全部 (224) 的特徵數量，再決定較佳之特徵數量，以下將說明五種方法如何排序特徵。

(一) 資訊獲利、獲利比率

這兩個方法顧名思義就是以資訊量中的資訊獲利來 評估屬性，資訊量定義為：「假設一個事件 S 有 n 種 結果，其發生機率分別為 P(v1), … P(vn)，則此事件發生後所得之資訊量定義如式 (1)：

2 1

( ) ( ) log ( )

n

i i

i

Entropy S P v P v

=

∑

− ⁽¹⁾

各種結果發生機率愈平均，所求資訊量也愈大，資訊量可以當作熵值 (entropy) 的指標，資訊量愈大，表示亂度愈大，以這樣的概念來解決特徵選取的問題。

資訊獲利是一種資訊理論的特徵選擇法，為集合 S 經 由屬性 A 分割後，所獲得的資訊量則為分割前的資訊 量減去分割後的資訊量，公式如式 (2)：

( , ) ( ) ^v ( _v)

v

Gain S A Entropy S S Entropy S

= −

∑

S ⁽²⁾

其中，|S| 代表集合 S 中所有資料的個數，|Sv| 為資料 集中特徵 A 等於 v 的樣本數，Gain(S, A) 即表示測試 前的資訊量減去測試後的資訊量，因為資訊愈小表示凌亂度愈小，所以要選擇測試後資訊量最小的特徵值，即選擇 Gain 最大的特徵值。獲利比率則是考慮當特徵值域太大時，雖然資訊獲利變大，但不一定是好的分割特徵，例如使用員工編號或身份證號碼當作特徵來分類時，其分割後的資訊量為 0，資訊獲利為最大。獲利比率定義為資訊獲利與分割資訊量 (split-info) 之比率，公式如式 (3)：

( , ) ( , )

Split-info( ) Gain S A Gain ratio S A

− = S (3)

其中 2

1

( ) - log

n

i i

i

S S

Split - info S

S S

=

∑

代表將集合 S 分成 n組 Si子集合[22, 23]。本研究針對每一欄位計算出該欄位進行分割所得到的資訊獲利及獲利比率進行排序，作為欄位選取之順序，以此方法從特徵集合中選取部份重要特徵。

(二) ReliefF

Relief演算法從訓練集中隨機選取 m 個樣本，透過所

選取的樣本與其同類和不同類的二個最近鄰居樣本的差異，求出每個樣本的各特徵與類別的相關性，再求平均值作為每個特徵的權重值，就得到每個特徵與類別的相關性。將 Relief 演算法擴展，不僅從同類與 不同類中各選出一個最近鄰居，而是選出 k 個最近鄰 居，求平均值得到每個特徵權重值，把特徵依權重值由大到小排序，即可依排序得到想要特徵數量的特徵，此為 ReliefF 演算法[24]。

ReliefF特徵選擇的輸入為選出的 m 個樣本及其所屬

類別，輸出量測特徵品質之權重向量 W[A]，A 為所有 特徵，其維度與特徵個數相同，表示特徵與類別的相關性。其演算法如下：

初始化：W[A]=0.0；

迴圈 i 從 1 至 m，執行

隨機選取 m 個樣本中的一個實例 Rⁱ；

找出與其同類(nearest hit Hj)的 k 個最近鄰居樣本；

對於每一類別 C ≠ class(Ri)使用群心法找出類別 C (nearest miss M_j(C))的 k 個最近鄰居樣本；對於每一 特徵 A，計算 W[A] (如式 4)：

1

1 ( )

( , , ) [ ] [ ]

( , , ( )) ( )

1 ( ( ))

i

k

i j

j

k

i j

C class R

diff A R H

W A W A

m k

diff A R M C P C

P class R m k

m k

=

≠

= − +

×

 

 − × 

 

×

∑

∑ ∑

(4)

取出下一個樣本，直到 m 個樣本都計算完成；

輸出 W[A]。

其中群心法為先計算某個群集之中心點，再找出非此 群集之 k 個最近鄰居。函式 diff(A, I1, I₂) = 0；當值 (I1, I₂) = 值 (I1, I₂)，否則 diff(A, I1, I₂) = 1。P(C) 與 P(class(R_i)) 分別為類別 C 與 class(Ri) 之事前機率，

求出各特徵與類別的相關性權重值 W[A]後，將其排 序，選擇所需的特徵數量。

(三) 支援向量機遞迴特徵消去法 (SVMRFE)

支援向量機 (SVM) 是貝爾實驗室研究人員 Vapnik [25]在對統計學習理論多年的研究基礎上發展起來的一種全新的機器學習演算法，SVM 是以 kernel 為基礎，有深厚的統計學習理論基礎，相較其它分類器像決策樹、類神經網路有些優點。大部份學習演算法是

(7)

採用經驗風險最小化，將訓練資料的誤差最小化。

SVM 期望能找出邊界 (margin) 最大的超平面 (hyper-plane)，將其方程式定義為 W^TX = −b，可寫為 W^TX + b = 0，其中 W^T＝{w1, w₂, …, w_n}為權重向量。

因此，兩支援超平面可寫為：

H1: w^Tx + b + δ^{, H2 : w}^T^{x + b}− δ

利用常數將 w、b、δ做 scaling，可將上兩式重寫為：

H1: w^Tx + b = 1, H2 : w^Tx + b = −1

從 H1 到原點的距離為 1−b/ w ，其中 |1−b|為 1-b 之絕對值、 w = w w⋅ 為 w 的歐基里德正規值 (Euclidean norm)。從 H2到原點的距離為 1− −b/ w ， 因此兩平面間的邊界為： 2/ w ，希望邊界最大，也 就是要最大化 2/ w ，即最小化 w / 2。利用以上的等式，在空間中的資料點必須滿足：

(w^T x_i) + b ≥ 1, if yi = 1 (w^T xi) + b ≤ −1, if yi = −1

可將以上二式改寫為 yi((w^T x_i) + b) ≥ 1，由以上說明，

可將 SVM 的問題歸納為：

.

min1 2

(( ) ) 1

T w b

T

i i

w w

y w x b



 + ≥



將此限制式解出即可得到我們所需要的超平面，過程中使用的 kernel 為 Gaussian kernel。支援向量機遞迴特徵消去法為 Guyon 等[26]於 2002 提出，其概念是先使用支援向量機演算法，計算在訓練模型中線性組合的特徵權重，再利用遞迴特徵消去法，將每個特徵權重平方後當作準則排序，刪除具最小準則之特徵，使用如此方法便可產生特徵子集合。

(四) OneR

OneR特徵選擇是以 OneR 分類器來評估特徵，它只用

一個特徵來分類，相當於一層的決策樹，其選擇特徵是以錯誤率來評估，選出分類錯誤率最低的特徵當作最重要的特徵，以此方式即可排序出特徵，其演算法如下[27]：

針對每一特徵 A

對於這個特徵之每一資料值 VA，產生如下規則：

計算每一特徵出現次數；

找出頻率最高之特徵 Ct；

當 A = VA，建立一規則，指派特徵值給 Ct；計算所有規則之錯誤率；

選擇具最小錯誤率之規則；

經過特徵排序後，接著執行特徵選擇，即從給定的特徵集合中，找到較佳之特徵子集合，由子集合產生的函數、訓練模型能夠達到充分的一致性[28]。Koller and

Sahami [29]認為給定某分類問題，一定可以找到特徵

子集合，使其產生的類別分佈無限接近訓練類別分

佈，該研究充分指出特徵選擇的可行性。評估特徵選擇效能最簡單的方式是以分類器的正確率為評估指標，一些研究宣稱假如目的是要最小化分類錯誤率，

而且所有特徵的錯誤成本相同時，分類器的預測正確率是最佳指標[30, 31]。將特定的學習演算法應用至特徵選擇程序中，然後依據特定搜尋演算法從原始訓練資料中挑出分類誤差最小之子集合，當作特徵選取結果，此種特徵選擇方式稱包裝式 (wrapper)；此外，

某些研究使用一些評分準則來排序特徵[32-34]，大部份是根據距離及資訊量來衡量，然後依排序結果將不重要特徵排除掉，此種特徵選擇方式稱過濾式

(filter)。雖然包裝式特徵選擇法可以由資料中挑選出

最佳的特徵子集合，但卻要耗費較高的計算成本，因此本研究使用過濾式特徵選擇法，以確保較佳的工作效率。

5. 分類

分類法是針對欲處理且未分類的資料集合，根據已知類別的物件集合，將欲處理資料依據其特徵去完成分類的過程，並冀望學習分類的規則，提供未來能自動分類之用；

亦即用「訓練資料」中的各種特徵，建構一個分類器，再使用訓練資料特徵相同、但資料內容不同的「測試資料」

驗證分類器之正確率。本研究在分類階段，評估過去辨識 ANA影像所使用的分類方法 (類神經網路、決策樹、樸素貝氏分類)。除此之外，由於支援向量機在許多領域分類具較佳效能，故本研究亦將支援向量機納入評估。

(一) 類神經網路

類神經網路的運作分為兩個階段來進行：(1) 訓練階段：用來調整網路內部各節點之間的鍵結值，以便輸入值經過網路計算之後能夠得到目標函數值；(2) 測試階段：主要用來驗證整個網路的精確性或者是用訓練完成的網路來進行事後的預測。類神經網路的內部架構分為三層，首先為輸入層，用來接收外來的資訊並將此資訊傳送到類神經網路中，以便進行運算，因此輸入值的類型必須要是數值的；第二為隱藏層，對於輸入層所接收到的資訊進行處理，不過使用者無法看見這個階段的處理流程；最後為輸出層，主要用來將隱藏層所處理過後的資訊，將之傳送至外界[22, 35, 36]。

類神經網路的訓練就是在調整鍵結值，使其變得更大或是更小，通常由隨機的方式產生介於+1 到−1 之間的初始值，鍵結值可視為一種加權效果，其值越大，

則代表連結的神經元更容易被激發，對類神經網路的影響也更大。

(二) 決策樹

決策樹是一種語意樹，與資料結構中的樹狀結構相仿，皆擁有根、節點及葉等結構，而每一節點有一分類條件，如「IF-THEN」的控制結構，利用測試結果

(8)

來決定資料將分類於此結點的哪一分支，並繼續作為分類的條件及最後的決策。Quinlan [37]提出 C4.5 演算法，其基本理論是改良他在 1979 年所提出的 ID3 演算法。C4.5 演算法是先建構一棵完整的決策樹，再針對每一個內部節點依使用者定義的錯誤預估率 (predicted error rate) 來作決策樹修剪的動作。另外，

C4.5演算法也和 ID3 相同，一樣是選擇資訊獲利最大的屬性為分割屬性。但是，兩者最大的差別在於 C4.5 特徵離散化的動作是由各內部節點動態決定的，亦即在不同的節點中，特徵值離散化的結果是不相同的。

建構決策樹的過程，以資訊獲利為依據，並選擇最大的資訊獲利值作為分類屬性。決策樹以此屬性的值分割成多個訓練子集合，而各子集合重複上述的步驟，

再從尚未被選為分類的屬性中找出最大的資訊獲利值，繼續作為分類的屬性，將決策樹分割成多個子樹直到不能再分類為止。

(三) 樸素貝氏分類

樸素貝氏分類主要是根據貝式定理，交換事前機率和事後機率，配合決定分類特性的各屬性彼此間是互相獨立條件機率的假設，來預測分類的結果。依數學學理來說明，即是根據訓練樣本，對於所給予測試訓練 資料中屬性值 (a1, a₂, a₃, …, a_n，假設一共有 n 個學習 概念) 的屬性 (A1, A₂, A₃, …, A_n)，a1_為A₁_{相對應的屬} 性值，指派具有最高機率值的類別 C (表示類別的集 合) 為目標結果，相關的演算法如下所述：

計算各個屬性的條件機率：

1 1 2 2

( _j| , , , _n _n)

P C=c A =a A =a A =a 貝氏定理：

1 2

1 2 1 1 1 2

( , , , | ) ( ) ( | , , , )

( , , , ) ( , , , | ) ( )

( , , , | ) ( ) ( , , , | ) ( )

n j j

j n

n

n j j

n n j j

P a a a c P c P c a a a

P a a a P a a a c P c

P a a a c P c P a a a c P c

=

= + +

假設屬性獨立，所以

1 2

1

( , , , | ) ( | )

n

n j i j

i

P a a a c P a c

=

∏

因此預測新訓練樣本所應歸屬的類別

1 2

1 1

arg max ( | , ,..., ) arg max ( ) ( | )

j

NB c C j n

n

c C j j

i

c P c a a a

P c P a c

=

= =

=

∏

其中 arg 為產生最大 P(cja1, a2, …, an) 的類別；ai為 表示與目標值相關聯的屬性值；cj為表示目標屬性；

P(cj) 為表示目標屬性值的機率；P(a1, a2, …, ancj) 為 表示每個屬性值在目標屬性 cj中出現的頻率。

根據訓練資料嘗試推估兩個機率：P(a1, a₂, …, a_ncj) 與 P(cj)，假設有 n 個與目標值相關聯的屬性值，那麼 就有 n(n-1) 個關聯性需要被考慮。假設與目標值相關 聯的屬性值，彼此為相互獨立，則從訓練樣本中選取 的目標值，即為目標值的機率 cNB，而 cNB即為個別屬性值機率的相乘積。

總而言之，只要樸素貝氏分類所及學習概念的屬性，

彼此間互相獨立的條件被滿足時，樸素貝氏分類所得 到的最大可能分類結果 cNB，與貝氏定理的最大可能 結果 P(CjA) 具有相同的意義。

6. 評估

在圖形辨識和機器學習的研究當中，通常會將資料分為訓練集和測試集，訓練集用來建立模型，測試集則為評估該模型，對於未知的樣本進行預測時的精準度；因此只有訓練集的資料才可以在模型的訓練過程中使用，然而測試集的資料必須要在模型完成後才拿來評估模型的好壞。

要將資料完整的分成訓練集和測試集，必須經由兩個準則：(1) 訓練集中的樣本數量必須充足，大致來說，至少大於總樣本數的 50%；(2) 兩組資料集合必須要從完整的集合中均勻取樣，主要的目的是希望能夠減少訓練集或測試集與完整資料集合之間的偏差[38]。

本研究分訓練集及測試集的方法是採用十折交叉驗證 (ten-fold cross validation)，即將資料分為十等分，9/10 做為訓練集，1/10 做為測試集，重複十次。訓練集之資料主要應用於圖 1 中監督式分類器之建構上，而測試資料主要是用來評估分類器之優劣，評估分類結果的指標有三：

(1) 正確率 (accuracy)；(2) Kappa 統計量；(3) 誤差均方根 (root of mean square error, RMSE)。正確率是用來計算分類系統對整體資料分類正確的比率，也就是被正確預測的資料個數佔全體資料個數的百分比，適用於資料的分佈平衡時，以表四說明，正確率公式如下：

a d accuracy

a b c d

= +

+ + +

當資料分佈不平衡時，則適合使用 Kappa 統計量，公 式如下：

1

o c

c

P P

Kappa P

= −

−

其中 Po為觀測一致性 (observed agreement)，表示前後測量結果一致的百分比，即實際為類別 A 預測也為類別 A 之機率，與實際為類別 B 預測也為類別 B 之機率和，計算公 式與正確率一致；Pc為期望一致性 (chance agreement)，表示前後測量結果預期相同的機率，預期機率為列總和乘以 行總和再除以總數，所以 Pc為實際為類別 A 預測也為類別 A 之預期機率與實際為類別 B 預測也為類別 B 之預期機率和：

(9)

表四類別預測結果表預測

類別 A 類別 B 總數類別 A a b a + b 實際類別 B c d c + d

總數 a + c b + d n = a + b + c + d

表五 Kappa統計量等級[39]

Kappa值吻合度

0.00~0.20 極低的吻合度 (slight)

0.21~0.40 一般的吻合度 (fair)

0.41~0.60 中等的吻合度 (moderate)

0.61~0.80 高度的吻合度 (substantial)

0.80~1.00 幾乎完全吻合 (almost perfect)

表六 ANA影像樣本分佈

ANA類別影像樣本數 ANA 類別影像樣本數

1. 廣泛型 479 4. 離散點點型 607

2. 大點點型 429 5. 周邊型 477

3. 小點點型 100 6. 核仁型 109

( ) ( ) / ( ) ( ) / /

Pc=a+ × +b a c n+ + × +c d b d n n

Kappa統計量是測量信度的指標，是在表現重複測量間的

一致性，Kappa 計算的結果介於 −1 到 1 之間，但通常落在 0~1，可分為五組來表示不同等級的吻合度，如表五所示。

第三個指標是 RMSE，是一種測量分類模型好壞之「誤 差」統計量，為觀察值與預測值之差值 (誤差) 平方的平均再開方根，此值將誇大離異值之存在，因此越小代表分類模型越好。

2 1/ 2

(O_i E_i)

RMSE n

 − 

=  



∑



其中 Oi為第 i 個觀察值，Ei為第 i 個預測值，n 為樣本數。

三、實驗結果與分析

本研究所使用之抗核抗體螢光顯影係由中部某教學

醫院提供，拍攝器材為 Olympus E-330，影像大小為 3162 像素 × 2352 像素，六類影像總計選取 2201 個細胞，其分佈情況如表六。影像之樣式判定是由三位專業技師與醫師共同認為該影像是屬於同一類，才將其歸類，若持不同意見之影像則不納入本研究。影像處理之程式以 Matlab

7.0.1開發，資料探勘軟體則是 Weka 3.4 [27]。本研究之影

像處理流程如圖 6，以下將逐一說明每個步驟：

實驗首先載入 ANA 影像，為計算影像中的灰階特徵，

需先將 RGB 影像轉為灰階影像 (圖 7)，轉換公式如下：

0.299 * 0.587 * 0.114 *

Gray= R+ G+ B

為使細胞邊緣檢測較為容易，在使用邊緣檢測器前，

會先對影像灰階值做對比度擴展，在此階段需要設定兩個 參數，門檻值 m 及控制轉換函數斜率的 E，本研究經實驗 將 E 設為 40，其對比有較佳結果，而門檻值 m 針對不同 影像，會設定不同數值，原則上影像較暗，m 值會較低，

反之亦然，到此為進行邊緣檢測之前處理。本研究評估的邊緣檢測方法有 Canny、小波轉換，Canny 檢測中高斯濾

(10)

54 技術學刊第二十八卷第一期民國一○二年

圖6 影像處理流程

圖7 RGB影像轉灰階影像

圖8 小波轉換後之影像

波器的標準差值設為 1，門檻值因子設 1.5，臨界值 T1、T2

是由系統自動計算提供 (T2為執行 Matlab 之 Canny 函數回 傳值乘以門檻值因子，T1 = 0.4*T2)。而小波轉換所採用的函數是 1910 年由 Haar 所提出的 Haar 函數[40]，Haar 函數不但結構簡單且運算快速，因此本研究採用 Haar 函數來對影像內容進行分析，經其轉換後之影像如圖 8。

使用邊緣檢測器得到的細胞輪廓常會有斷點、不完整的情況，在此使用形態分析來改善，其運算步驟先將輪廓影像膨脹 (圖 9(a))，再將影像中所有的區域及洞填滿 (圖

9(b))。由於在影像邊緣的細胞是不完整的，所以將其去除

(圖 9(c))，去除後再針對餘下的細胞進行侵蝕 (圖 9(d))，即完成形態分析的步驟。為了比較邊緣檢測的結果與原圖的吻合度，將檢測出之邊緣及原圖重疊，以利觀察 (圖 10)。

圖9 形態分析

圖10 檢測之邊緣與原圖對照圖

經過形態分析改善邊緣後，自行選取有興趣之區域 (region of interest, ROI) 得到一完整細胞，即可透過紋理分析結合灰階共生矩陣 (gray level co-occurrence matrix)、遊程矩陣 (run-length matrix) 分別計算出 180 個與 44 個，合計 224 個灰階紋理特徵。由於各特徵之值域差異甚大，為了避免值域較大特徵會減弱值域較小特徵的影響，故在進行下一階段之前，本研究針對每個特徵做標準化，使所有特徵的值介於 −1 到 1 之間。由於目前特徵維度偏高，將影響機器學習速度，且特徵集合內可能含有無關或相似的特