基于多模態特征的醫學圖像聚類方法*

王保加，潘海為，謝曉芹，張志強，馮曉寧

哈爾濱工程大學 計算機科學與技術學院，哈爾濱 150001

1 引言

隨著生活水平的提高，醫療健康問題受到人們越來越多的關注，隨之帶來了醫療數據的飛速增長。醫療大數據可以有效地應用在疾病的預測分析、規范分析和基因組學研究方面[1]。醫學圖像又是醫療大數據的重要組成部分，腦部CT圖像由于在檢測腦梗塞、腫瘤、鈣化、出血和骨創傷時效果顯著，從而被廣泛使用[2]，而圖像聚類方法在圖像早期處理和深度挖掘中發揮著重要作用。本文旨在分析腦CT圖像的特征，實現對腦CT圖像的有效聚類，以輔助醫生診斷并挖掘出潛在的信息。

現有的圖像聚類方法主要有兩種：一種是基于文本的圖像聚類；另一種是基于內容的圖像聚類。基于文本的圖像聚類方法需要人工標注特征描述[3]。基于內容的圖像聚類方法則通過提取圖像特征，然后計算特征之間的相似度來實現聚類[4]，因此特征的選取和相似性的度量決定了聚類結果的準確性。按照特征的表現方式，圖像特征可分為顏色特征、形狀特征、局部特征和紋理特征等。顏色特征如顏色直方圖（color histogram）[5]和顏色矩（color moment）[6]常應用于彩色圖像、遙感圖像等多光譜圖像。圖形的骨架結構[7]和邊緣檢測（edge detection）[8]能夠有效地表達簡單圖像的形狀特征。局部特征如尺度不變特征轉換（scale-invariant feature transform，SIFT）[9]可以解決兩幅圖像之間發生平移、旋轉、視角變換、光照變換等情況下的匹配問題。紋理特征能精簡有效地描述圖像主體的信息，對相似結構圖像的特征提取可以取得良好效果。常用的紋理特征有局部二值模式（local binary pattern，LBP）[10]、灰度共生矩陣（graylevel co-occurrence matrix，GLCM）[11]、Gabor小波[12]、方向梯度直方圖（histogram of oriented gradient，HOG）[13]等，紋理特征在圖像聚類[14]、人臉識別[15]、圖像檢索[16]等方面應用廣泛。

然而醫學圖像特征提取和普通圖像相比有很大的難點，醫學圖像的相似性度量涉及到病理學和CT成像相關知識，導致普通圖像的特征提取方法無法完全表達醫學圖像的全部信息。比如對于在不同區域發生相同病變的醫學圖像，它們在病理學上應歸屬同類別的病變圖像，但現有的特征提取方法忽略了這種高層語義信息，不僅不能突顯不同位置尤其是病變區域的重要性，還會因為正常區域和病變區域的相互干擾造成圖像相似性比較不準確。同時醫生對CT圖像診斷時，也更加關注具有診斷意義的ROI（region of interest）區域，而忽略細微紋理的差異。因此提取出能夠反映人類視覺認知的醫學ROI特征，在醫學圖像的聚類中具有重要的意義。

基于圖像聚類的研究已經取得了很多進展。如Zhan等人提出了一種基于圖熵的醫學聚類算法[17]；Zhu等人使用了顏色、紋理、邊緣的全局和局部特征進行圖像匹配[18]；Pan等人提出提取圖像中的ROI區域來表征一幅圖像，但其沒有考慮圖像中不同ROI區域的重要程度的差異[19]；Mour?o等人融合文本和圖像數據實現了對醫學圖像的聚類和檢索[20]；Liang等人分別把顏色、HOG和SIFT特征向量拼接成新向量完成圖像聚類[21]。這些圖像聚類方法都只是對圖像的一類特征進行聚類，不能對圖像的多模態特征選取不同特征的權值，進而達到特征融合的目的。

針對以上研究存在的問題，本文提出了一種融合圖像紋理特征和醫學圖像特有領域知識的形態學特征的多模態特征抽取方法。首先利用紋理特征在相似紋理圖像上的良好區分效果，進而融入醫學圖像領域知識，結合醫學圖像ROI在病變圖像的區分能力，實現了對復雜腦部CT圖像集的聚類，并使用了多核譜聚類算法驗證了聚類的效果。

本文組織結構如下：第2章介紹多模態特征醫學圖像聚類的總體框架；第3章研究紋理特征和形態學特征提取方法；第4章討論基于多模態特征的聚類方法；第5章展示實驗結果和相關分析；最后總結全文。

2 多模態特征聚類總體框架

本文根據病理學和CT成像相關知識，提出了提取醫學圖像中具有診斷意義的ROI區域的形態學描述表達圖像。CT圖像是以不同的灰度值來反映不同密度的器官組織對X光射線的吸收程度[22]。低密度區域的組織對X光射線吸收少，在CT圖像中顯示為黑影；高密度區域的組織則顯示為白影。圖1（a）為正常人的腦CT圖像；圖1（b）為腦出血CT圖像，白影為腦出血位置；圖1（c）為腦梗塞CT圖像，圖像左邊的黑影和中間腦室的深色黑影為腦梗塞位置，腦梗塞通常是黑影區。由于梗塞病變的疊加作用，圖1（c）中腦室位置成像比正常人腦室灰度值更低。醫生在對CT圖像診斷時，更加關注這些具有診斷意義的黑影和白影區域。因此本文使用具有醫學領域知識的黑影和白影區域作為ROI區域。

Fig.1 Different brain CT images圖1 不同腦CT圖像

本文分別從紋理角度和ROI形態學角度提取醫學圖像集的特征，圖2為多模態特征聚類方法總體流程。一方面對于紋理特征，分別選取了基于局部紋理特征的LBP、基于全局紋理特征的GLCM和基于形狀紋理特征的邊緣方向直方圖（edge direction histogram，EDH）3種特征來表達醫學圖像全局的細粒度底層語義信息，通過本文提出的基于紋理特征的提取方法（feature extraction based on texture，FET）構造紋理的特征矩陣MatrixT。另一方面本文提取出醫學圖像中具有診斷意義的ROI區域，從形態學角度分析ROI的灰度級、大小、位置、方向、形狀等來表達醫學圖像局部的粗粒度高層語義信息。首先通過ROI區域提取（ROI extraction based on gray and area，REGA）方法提取到ROI，然后使用圖像ROI向量化（ROI vectorization，ROIV）方法實現圖像同維向量表達，最后得到圖像集基于形態學的特征矩陣MatrixR。根據兩種特征的互補作用關系，本文定義了一種多模態相似性函數（multimodal similarity function，MSF）計算兩種特征的相似度矩陣，并對兩類特征通過多核譜聚類（multi-kernel spectral clustering，MKSC）算法對圖像聚類，并得到兩種特征的特征融合權值。

Fig.2 Flow chart of multimodal features clustering圖2 多模態特征聚類總體流程圖

3 特征提取

3.1 紋理特征

本文考慮到不同的紋理特征的適用性，由于CT圖像對光照極其敏感，使得同一個人在不同醫院不同設備下的CT圖像存在差異，而LBP特征具有旋轉不變性和灰度不變性優點，因此使用LBP特征可以有效地降低光照的影響；GLCM使用統計方法表達整幅圖像，可以反映圖像全局的紋理粗細度及方向性對比變化；同時腦部各切片的CT圖像由于掃描的截面不同，造成相鄰切片CT圖像形狀相似，但大小不一，EDH可以通過檢測紋理的邊緣方向提高各切片之間的關聯度。據此本文選取這3種紋理特征描述CT圖像，以下將對3種特征進行介紹，并提出基于紋理的醫學圖像特征提取方法。

3.1.1 紋理特征描述

（1）LBP

LBP特征通過比較中心像素點和鄰域半徑內像素點的灰度值，計算出中心像素點二進制編碼序列，再統計不同LBP序列的數目來表達圖像的LBP紋理特征FL。對于一個像素點的LBP值定義為：

其中，LBPP,R表示在中心像素點半徑R上的P個采樣點；gc為鄰域中心像素的灰度值；gi為半徑范圍內第i個像素的灰度值。

例1如圖3所示為LBP值計算過程：實驗中選擇半徑1上的8個像素LBP8,1進行采樣。即每個采樣窗口為3×3的方陣，左邊方陣中心像素灰度值gc為88，gc半徑1內的8個采樣點gi分別為10、60、189、55、211，75、58、96（從方陣左上角順時針遍歷），計算S(gi-gc)得到右邊方陣8位0和1組成的LBP序列：（10010100）2=148，148即為中心點的LBP值，總共形成最多可以生成28=256種不同狀態的LBP值。

Fig.3 Calculate LBP圖3 LBP值計算

統計一張圖像中不同LBP值的數目作為圖像的特征向量，則一張圖像可由256維向量表示。Ojala等人提出了均勻LBP模式，均勻LBP定義為0-1跳變不超過2次的編碼。研究者發現圖像中均勻LBP在圖像所有LBP模式中占比達90%以上，因此提取58類均勻LBP和非均勻模式組合成59維向量進行降維處理[23]。本文選用均勻LBP提取圖像特征。

（2）GLCM

GLCM是像素距離和角度的矩陣。本文使用灰度共生矩陣的角二階矩、對比度、相關性和熵值來提取GLCM的紋理特征

角二階矩（angular second moment，ASM）是灰度共生矩陣元素值的平方和，反映了圖像灰度分布均勻程度和紋理粗細度。

對比度（Contrast）反映了圖像的清晰度和紋理溝紋深淺的程度，紋理溝紋越深，其對比度越大。

相關性（Correlation）度量了空間灰度共生矩元素在行或列方向上的相似程度。

熵（Entropy）是圖像所具有的紋理信息量的度量，若圖像沒有任何紋理，則熵值接近為0。

其中，Pij表示在θ方向上，相隔距離d的一對像素分別具有灰度值i和j出現的概率；L表示灰度級。

例2分別從 0°、45°、90°、135°共4個方向計算出圖4（a）的GLCM的16個相關值，ASM=[0.874,0.848,0.877,0.849]，Contrast=[0.460,0.655,0.423,0.626]，Correlation=[0.953,0.932,0.957,0.936]，Entropy=[0.204,0.191,0.205,0.192]。

（3）EDH

EDH提取了圖像形狀邊緣紋理特征。如圖4是EDH的提取過程。首先將圖像4（a）進行canny邊緣檢測，得到邊緣紋理圖像4（b）；然后對邊緣紋理上的像素點I分別沿x軸、y軸計算邊緣方向梯度值dx和dy，則像素點I的邊緣方向如圖 4（c）對 [-180°,180°]角度區間每隔 10°劃分一組，一共分成36組邊緣方向區間，統計像素點落在每組區間的個數，最終生成36維邊緣方向向量FE。

Fig.4 EDH extraction圖4EDH提取圖

3.1.2 FET方法

圖像整體地提取特征容易忽略局部重要的紋理信息，因而對圖像劃分子圖分別提取特征被廣泛使用[16，25]。FET算法的基本思想是把每一張圖像分割成大小相同的區域，首先對每個區域分別求取3種紋理特征向量，再對所有區域的特征向量進行合成，最后使用主成分分析（principal component analysis，PCA）進行降維，算法輸出紋理特征矩陣MatrixT。

對圖像集IM={IM1,IM2,…,IMN}，第i張圖像IMi分割得到r個大小為s×t的矩形分割區域。對于第j個矩形分割區域分別使用3種特征共111維向量表示此區域則圖像IMi紋理特征向量表達為表示為D(111×r)行1列，每張圖像通過PCA降到Z維，最后輸出N×Z的紋理特征矩陣MatrixT。具體過程如算法1所示。

算法1FET算法

輸入：圖像集IM={IM1,IM2,…,IMN}。

輸出：圖像集IM的紋理特征矩陣MatrixTN×Z。

1.ForIM的每一張圖像IMi

2.ForIMi的每一個分割區域j

3.計算第j區域的添加到向量數組

4.End for

6.End for

7.ForXD×Z的每一列Xi

8.計算Xi的平均值E(Xi)；

9.End for

10.計算協方差矩陣：

11.求出協方差矩陣Z個最大的特征值及對應的特征向量UZ×D=[μ1,μ2,…,μZ]T；

12.返回MatrixTN×Z=(UZ×D?XD×N)T；

算法1中第1～7步求取每張圖像特征向量XD×N，時間復雜度為O(N)，第8～12步計算協方差矩陣的特征值時間復雜度為O(N3)，算法總的時間復雜度為O(N3)。

3.2 形態學特征

醫學圖像ROI的提取是為了獲取醫學圖像中具有診斷意義的像素區域。本節給出形態特征ROI提取的方法、ROI相關描述和圖像ROI向量化的方法。

3.2.1 REGA方法

本文使用具有醫學領域知識的黑影和白影區域作為ROI區域。以下是ROI提取的過程，定義1給出了一個圖像集中ROI的相關概念。

定義1給定數量為N的圖像集IM={IM1,IM2,…,IMN}，設第i張圖像IMi擁有Ki個ROI，則圖像集第i張圖像的第j個ROI定義為Rij，{Rij∈IMi|1≤i≤N,0≤j≤Ki}。

本文提出了ROI區域提取算法REGA，算法的基本思想是根據黑影灰度閥值BGrayε、白影灰度閥值WGrayε和面積閥值Areaε的限制分割出最佳的ROI區域，具體過程如算法2所示。

算法2REGA算法

輸入：圖像集IM中的第i張圖像IMi；灰度閥值BGrayε和WGrayε；面積閥值Areaε。

輸出：圖像的ROI集合{Ri1,Ri2,…,RiKi}。

1.IMi使用sobel算子進行邊緣檢測、高斯濾波噪聲抑制、形態學開閉操作平滑去噪；

2.IMi提取區域邊緣并進行二值化處理，得到分割圖像IMB；

3.For圖像IMB分割區域

4.IfIMB分割區域面積＜Areaε

5. ThenIMB分割區域值設置為0；

6.End if

7.Else ifIMi邊緣灰度值＜BGrayε或邊緣灰度值＞WGrayε

8. ThenIMB分割區域值設置為1；

9. ElseIMB分割區域值設置為0；

10. End if

11.End for；

12.標記IMB為1的連通區域，返回圖像IMi的ROI集合 {Ri1,Ri2,…,RiKi}；

算法2中輸入圖像IMi通過第1～2步把不同灰度值像素分割成多個區域的二值化圖像IMB；第4～5步設置了面積閥值Areaε過濾掉部分小的紋理區域；第6～8步設置灰度閥值BGrayε和WGrayε，分別用于提取黑影和白影區域，通過調整3個值來獲得最佳ROI區域；最終得到具有Ki個連通區域的IMB圖像，即為所要提取的ROI。算法2的時間消耗在對圖像像素點的遍歷上，因此對圖像集算法2總的時間復雜度為O(N)。

圖5為ROI區域提取圖。圖5（a）為腦出血圖像；圖5（b）為腦出血圖像提取的8個黑影ROI；圖5（c）為腦出血圖像提取的4個白影ROI；圖5（d）紅色邊緣為原圖上重新標注的黑影和白影位置，這些ROI很好地反映了圖像整體的病理學特征。

對于提取到的ROI，從形態學的特征上可以很好地描述不同ROI的特征。本文使用ROI的面積、ROI的位置、ROI的離心率、ROI的周徑比、ROI的圓形度、ROI的方向和ROI黑影白影標記來表達一幅圖像。對于一個ROI區域R，具體定義如下：

R的面積（Area）表示為AreaR=M，其中M為ROI區域的像素點的總個數。

R的重心位置（Location）表示為LocR=(Centroidx,Centroidy)，其中Centroidx、Centroidy表示該區域的二維重心坐標值。

Fig.5 ROI area extraction圖5 ROI區域提取圖

R的離心率（Eccentricity）表示為其中LAL表示與R具有相同標準二階中心矩的橢圓的長軸長度（long axis length，LAL），SAL表示與R具有相同標準二階中心矩的橢圓的短軸長度（short axis length，SAL）。

R的周徑比（perimeter radius ratio，PRR）表示為其中Per表示R的周長，Rad表示與R具有相同面積的圓的半徑。

R的圓形度（Cularity）表示為其中Ave表示R的重心到邊界各點的平均距離，σ表示R的重心到邊界各點距離的方差。

R的方向（Orientation）表示為OriR=θ，其中θ為與R具有相同標準二階中心矩的橢圓的長軸與X軸正方向的夾角。

R的灰度標記（gray label，GL）在病理學上是兩種不同的標記，這兩種R的相似度為0：

因此任意的兩個ROI區域Ri可由一個k=7維向量和灰度標記表示。對于任意的兩個ROI區域Ri和Rj，距離計算為：

3.2.2 ROIV方法

ROI圖像向量化的目的是把由多個ROI表示的圖像轉化成數值向量的表現形式，本文提出ROIV算法的基本思想是使用詞袋模型（bag of words，BoW）對擁有不同ROI數目的圖像實現同維向量表示，ROI詞袋模型的描述如定義2所示。本文考慮到不同ROI的重要程度差異，使用具有顯著性表征ROI重要性的指標，即ROI的灰度級和ROI的面積對ROI單詞賦予不同的權重。顯然ROI的灰度級和正常灰度級相差越大表示此ROI越重要（越可能是病變位置），ROI的面積越大表示此ROI也越重要。據此本文提出了（ROI term weight-inverse document frequency，RTWIDF）方法用于對詞典生成的數值向量特征加權，算法輸出圖像集的形態學特征矩陣MatrixR。

定義2設一個ROI單詞（ROI term，RT）是一個表示這個ROI區域的7維向量，則對于IM圖像集，ROI詞典（ROI document，RD）為ROI單詞的集合，RD={R11,R12,…,R1K1,…,RN1,RN2,…,RNKN}，其中N表示圖像的數目。IMi的RD可由Ki個RT表示。圖像集的總數（sum of RD，SR）為

對于圖像集的RD，利用式（7）計算距離并使用K-Means算法聚成C類，則圖像IMi可由向量Vi=(ni1,ni2,…,nic)表示，nij表示聚類中IMi在第j類中RT數目。

定義3R灰度級（gray level，GL）為R內所有像素點灰度的平均值，其中M為R中像素點的總數，Ii表示第i個像素點的灰度值。

定義4圖像集IM的正常灰度級（normal gray Level，NGL）為正常圖像樣本集ROI灰度級的平均值，，其中SR為正常圖像提取到的ROI數目，GLi為第i個ROI區域Ri的灰度級。

定義5對于一個ROI區域R，R的灰度級權重（gray level weight，GLW）為：

定義6對于圖像集IM中圖像IMi的第j個ROI區域Rij，Rij面積權重（area weight，AW）為：

RTW-IDF加權：

RTW：對于圖像集IM中圖像IMi的第j個ROI區域Rij，Rij權重（ROI weight，RW）為：

其中，ρ為大于1的常數，ROI灰度的影響要比面積更顯著。

圖像IMi的向量Vi的第j維中包含的單詞權值（term weight，TW）和為為聚類中IMi在第j類中RT數目，RWj為第j個RT權值。

IDF：圖像IMi的向量Vi的第j維中逆文獻頻率（inverse document frequency，IDF）為，其中N為圖像集中圖像的數量，nj是第j類RT包含的圖像數目。

最終圖像IMi特征向量權值化后的C維特征向量表示為VWi=(VWi1,VWi2,…,VWiC)，其中VWij=RTWij?IDFij。具體的過程如算法3所示。

算法3ROIV算法

輸入：圖像集IM中的ROI集合RD={R11,R12,…,R1K1,…,RN1,RN2,…,RNKN}聚類中心數C。

輸出：圖像集IM的形態學特征矩陣。

1.初始化N行C列的全零矩陣MatrixRN×C；

2.利用式（7）對RD進行K-Means聚成C類；

3.由第1步計算每一個類的聚類中心向量特征向量U={u1,u2,…,uc}，單詞數目n1,n2,…,nc，并對單詞數目使用IDF的定義計算得到idf數組；

4.ForIMi的每一張圖像

5.ForIMi的每一個ROI區域Rij

6. 歐氏距離計算Rij最近U的聚類中心，標記為類別p添加數組q，式（9）計算Rij在圖像IMi中的ROI權值rw；

7. 計算rw*idf[p]添加到數組w作為Rij在特征矩陣第p維的分量；

8.End for

9.For數組q的每一個元素q[i]

10.取回類別標記p=q[i]；

計算MatrixRip+=w[i]；

11.End for

12.End for

13.返回MatrixRN×C；

算法3中第2～3步使用K-Means對ROI單詞進行聚類，時間復雜度為O(|RD|×C×T)，|RD|為ROI單詞的個數，T為迭代次數；第4～15步是對圖像向量化的過程，其中需要遍歷整個圖像集和每張圖像的所有ROI，時間復雜度O(N×K)，K為一張圖像的ROI上限，由于C×T大于K，從而對圖像集算法3總的時間復雜度為O(|RD|×C×T)。

4 基于多模態特征的聚類算法

4.1 MSF方法

相似性度量選取對聚類結果有很大的影響，高斯核函數常被用于度量兩張圖像的相似度：

其中，dist(xi,xj)表示數據xi和xj歐氏距離；σ為尺度參數。為自適應確定σ，文獻[26]定義了相似度函數：

其中，σi表示數據xi與第L(L=7)個近鄰點的距離。

第3章利用特征提取算法得到了圖像集IM的紋理特征向量矩陣MatrixT和形態學特征向量矩陣MatrixB。考慮到兩種特征在圖像相似度上的互補作用，本文重新定義了多模態相似性度量函數MSF：

其中，φi表示圖像IMi分別使用紋理特征和形態學特征前L個近鄰點集合中公共交集的圖像數目。設A表示IMi在紋理特征度量下前L個近鄰圖像的集合，B表示IMi在形態學特征度量下前L個近鄰圖像的集合，則兩個集合中相同的圖片數量φi=A?B。

通過本文定義的MSF計算紋理特征向量矩陣MatrixT得到N×N的紋理特征相似度對稱矩陣同理計算形態學特征向量矩陣MatrixR得到形態學特征相似度對稱矩陣SMatrixR。

4.2 MKSC方法

多核學習（multi-kernel learning，MKL）是由Lanckriet等人提出的[27]，在此基礎上Huang等人又把多核學習方法擴展到非監督的聚類領域[28]，多核學習可以自適應地完成對不同特征分配權值。

由圖的譜理論[29]可知，對于給定的N個數據集X={x1,x2,…,xN}，當找到滿足式（12）的向量集U={u1,u2,…,uK},ui∈RK時，數據集X聚成K類。

其中，Wij表示數據xi和xj的相似度。

對于4.1節中數據集兩種特征相似度矩陣，給定紋理特征權值為α1，形態學特征權值為α2，且α1+α2=1，則式（12）轉變為求解帶約束條件的式（13）：

MKSC算法的基本思想是：求解式（13）時，首先初始化α1=α2=0.5；然后迭代以下過程，固定α1和α2，式（13）轉化為關于α1、α2的帶線性約束的二次規劃問題；固定參數α1和α2，式（13）轉化為關于α1和α2的帶線性約束的非線性最小化問題。交替求解上述兩個問題，直到滿足收斂條件，當收斂時結束迭代，求解到α1和α2并輸出聚類結果。具體過程如算法4所示。

算法4MKSC算法

輸入：圖像集IM={IM1,IM2,…,IMN}，相似度矩陣SMatrixT、SMatrixR，聚類種類K。

輸出：K個聚類類別 {Y1,Y2,…YK}，特征權值α1、α2。

1.初始化特征權值α1=α2=0.5；

2.Repeat

3.構造相似性度量矩陣：

5.計算L最小的K個特征值λ1,λ2,…,λK所對應的特征列向量u1,u2,…,uK；

7.更新α1、α2；

8.Until收斂

9.構造矩陣U={u1,u2,…,uK}，ui∈RK得到矩陣V，其中

10.對V行向量使用K-Means聚類算法聚成K類，得到聚類結果{Y1,Y2,…,YK}，如果V的第i行屬于j類，則IMi屬于Yj；

11.返回{Y1,Y2,…,YK}，α1，α2；

算法4中主要包括第2～8步迭代求解特征權值和第9～10步K-Means兩個過程，前一個過程的時間復雜度為O(N3T1)，T1為迭代求解次數，后一個過程的時間復雜度為O(NKT2)，T2為迭代次數，因此總的時間復雜度為O(N3T1)。

5 實驗

5.1 數據集和環境

實驗中使用1 600張某醫院的真實腦部CT圖像，圖像包括正常圖像、腦出血圖像、腦梗塞圖像。硬件環境：AMD A4-3300M@1.9 GHz，6.0 GB RAM，Windows 7系統。軟件環境：Matlab 2014a。

5.2 評價標準

本文采用標準化互信息（normalized mutual information，NMI）作為評價標準[30]。對于具有K個真實分類集T和聚類結果集C的NMI：

其中，I(T；C)為互信息；H為信息熵：

其中，n為樣本總數目；N(ti)為分類集T中第i分類的數目；N(cj)為分類集C中第j分類的數目；N(ti,cj)為在T中第i分類和C中第j分類同時出現的數據數目。

5.3 實驗設置和結果

紋理特征提取時需要對圖像進行分割，不同的分割方式會有不同的結果，通常分割后圖像的總數小于10[25]。基于此本實驗比較了5種不同的圖像分塊分割方式：1×1（不分區），2×2，3×3，2×4和4×2。表1為不同分塊模式下的NMI值，當分割為4×2的模式時，能取得最好的效果，因此本實驗使用4×2的模式進行圖像分塊。

Table 1 NMI values of different segmentation methods表1 不同圖像分塊模式的NMI值

本文還需要對圖像ROI向量化時向量維數C值進行設置，太小的C值會把不同類別的ROI聚成一類，太大的C值會使聚類太分散，造成圖像區分不明顯。圖6為不同C值下使用形態學特征圖像聚類的NMI值。當C取值在28～31附近時NMI值趨于最大且穩定于0.58，因此本文以下實驗C都設置為30。

Fig.6 NMI in different C values圖6 不同C值的NMI

圖7為正常圖像（normal）和病變圖像（abnormal）不同比例normal/abnormal下紋理特征的NMI值。當比值為0時表示圖像全部為病變圖像，當比值為1時表示全部為正常圖像。從實驗結果中可以看出，NMI值隨著正常圖像比例的增加而增大，原因是紋理特征受病變圖像的干擾會導致聚類的準確率降低；同時由4種曲線可以看出，3種紋理特征融合后的多模態特征可以取得比單個特征更好的效果。

Fig.7 Comparison of different texture features圖7 不同紋理特征比較

本文對比了4種圖像聚類方法，包括SIFT、SOM（self-organizing maps）、Zhan（2016）[14]以 及 Liang（2016）[21]。如圖8所示，本文首先使用了傳統的SIFT特征點詞袋模型和基于圖像像素聚類的SOM方法，SIFT特征點聚類方法隨著圖像的增加NMI值增大，這是由于圖像數目的增加使得圖像詞典更加完善，SOM聚類的NMI值在0.40和0.45之間波動，但兩種方法的NMI值都在較低的水平，可以看出這兩種傳統的方法不能很好地應用在醫學圖像的聚類上；Zhan（2016）使用圖熵聚類，在低數量級能取得比本實驗更好的效果，但是隨著數量的增加實驗NMI值降低；Liang（2016）使用了半監督的譜聚類進行圖像聚類，它通過訓練構造一個新的相似性度量函數，實驗對訓練的結果依賴度高。

Fig.8 Comparison of different clustering methods圖8 不同圖像聚類方法對比

本文提出的醫學圖像聚類方法不需要復雜的訓練過程，實驗的NMI值雖然在低數量級下低于Zhan（2016）的方法，但隨著圖像數目的增加實驗結果呈現增長趨勢。這是由于圖像集越大，醫學圖像的ROI類型越多越全面，對醫學圖像的特征描述更準確。本實驗的方法整體優于Liang（2016）方法。

本文提出的多模態特征的醫學圖像聚類方法的部分聚類結果如圖9所示，其中圖9（a）為正常圖像，圖9（b）為腦出血圖像，圖9（c）為腦梗塞圖像。

6 結束語

本文提出了一種融合圖像紋理特征和醫學圖像特有領域知識的形態學特征的多模態特征抽取方法，旨在解決傳統特征提取方法不能有效地表達醫學圖像特征導致聚類不準確的難題。本文提出了醫學圖像的ROI提取算法和ROI形態學特征描述，使用特征融合的多核譜聚類算法完成圖像聚類。實驗結果表明此方法可以減小非正常醫學圖像和正常醫學圖像聚類的相互干擾，實現復雜醫學圖像聚類。

Fig.9 Result of MKSC image clustering圖9 MKSC方法圖像聚類結果

[1]Adamson D.Big data in healthcare made simple:where it stands today and where it's going[EB/OL].(2016)[2017-11-22].https://www.healthcatalyst.com/big-data-in-healthcaremade-simple.

[2]Galloway R L.Introduction and historical perspectives on image-guided surgery[M]//Golby A J.Image-Guided Neurosurgery.Amsterdam:Elsevier North-Holland,Inc,2015:1-22.

[3]Demerdash O E,Kosseim L,Bergler S.Text-based clustering of the ImageCLEFphoto collection for augmenting the retrieved results[C]//LNCS 5152:Proceedings of the 8th Workshop of the Cross-Language Evaluation Forum Advances in Multilingual and Multimodal Information Retrieval,Budapest,Sep 19-21,2007.Berlin,Heidelberg:Springer,2008:562-568.

[4]Yang Lei,Tian Jun,Wu Dapeng.Content-based image authentication by feature point clustering and matching[J].Security and Communication Networks,2012,5(6):636-647.

[5]Swain M J,Ballard D H.Color indexing[J].International Journal of Computer Vision,1991,7(1):11-32.

[6]Stricker A M,Orengo M.Similarity of color images[C]//SPIE 2420:Proceedings of the Storage and Retrieval for Image and Video Databases III,San Diego/La Jolla,Feb 5-10,1995:381-392.

[7]Shen Wei,Wang Yan,Bai Xiang,et al.Shape clustering:common structure discovery[J].Pattern Recognition,2013,46(2):539-550.

[8]Zitnick C L,Dollár P.Edge boxes:locating object proposals from edges[C]//LNCS 8693:Proceedings of the 13th European Conference on Computer Vision,Zurich,Sep 6-12,2014.Berlin,Heidelberg:Springer,2014:391-405.

[9]Lowe D G.Distinctive image features from scale-invariant keypoints[J].International Journal of Computer Vision,2004,60(2):91-110.

[10]Ojala T,Pietik?inen M,Harwood D.A comparative study of texture measures with classification based on feature distributions[J].Pattern Recognition,1996,29(1):51-59.

[11]Haralick R M,Shanmugam K S,Dinstein I.Textural features for image classification[J].IEEE Transactions on Systems,Man,and Cybernetics,1973,3(6):610-621.

[12]Daugman J G.Uncertainty relation for resolution in space,spatial frequency,and orientation optimized by two-dimensional visual cortical filters[J].Journal of the Optical Society of America A Optics&Image Science,1985,2(7):1160-1169.

[13]Dalal N,Triggs B.Histograms of oriented gradients for human detection[C]//Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition,San Diego,Jun 20-25,2005.Washington:IEEE Computer Society,2005:886-893.

[14]García-Villalba L G,Orozco A L S,Corripio J R.Smartphone image clustering[J].Expert Systems with Applications,2015,42(4):1927-1940.

[15]Tan Xiaoyang,Triggs B.Enhanced local texture feature sets for face recognition under difficult lighting conditions[J].IEEE Transactions on Image Processing,2010,19(6):1635-1650.

[16]Subrahmanyam M,Maheshwari R P,Balasubramanian R.Local tetra patterns:a new feature descriptor for contentbased image retrieval[J].IEEE Transactions on Image Processing,2012,21(5):2874-2886.

[17]Zhan Yu,Pan Haiwei,Han Qilong,et al.Medical image clustering algorithm based on graph entropy[C]//Proceedings of the 12th International Conference on Fuzzy Systems and Knowledge Discovery,Zhangjiajie,Aug 15-17,2015.Piscataway:IEEE,2015:1151-1157.

[18]Zhu Jianke,Hoi S C H,Lyu M R,et al.Near-duplicate keyframe retrieval by nonrigid image matching[C]//Proceedings of the 16th International Conference on Multimedia,Vancouver,Oct 26-31,2008.New York:ACM,2008:41-50.

[19]Pan Haiwei,Li Jianzhong,Zhang Wei.Incorporating domain knowledge into medical image clustering[J].Applied Mathematics and Computation,2007,185(2):844-856.

[20]Mour?o A,Martins F,Magalh?es J.Multimodal medical information retrieval with unsupervised rank fusion[J].Computerized Medical Imaging&Graphics,2015,39:35-45.

[21]Liang Jianqing,Han Yahong,Hu Qinghua.Semi-supervised image clustering with multi-modal information[J].Multimedia Systems,2016,22:149-160.

[22]Kitamoto A.Data mining for typhoon image collection[C]//Proceedings of the 2nd International Workshop on Multimedia Data Mining,San Francisco,Aug 26,2001:68-77.

[23]Ojala T,Pietik?inen M,M?enp?? T.Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2002,24(7):971-987.

[24]Lu C S,Chung P C,Chen C F.Unsupervised texture segmentation via wavelet transform[J].Pattern Recognition,1997,30(5):729-742.

[25]Chen Yarui,Yang Jucheng,Wang Chao,et al.Multimodal biometrics recognition based on local fusion visual features and variational Bayesian extreme learning machine[J].Expert Systems withApplications,2016,64(C):93-103.

[26]Zelnik-Manor L,Perona P.Self-tuning spectral clustering[C]//Proceedings of the 2004 Conference on Neural Information Processing Systems,Vancouver,Dec 13-18,2004:1601-1608.

[27]Lanckriet G R G,Cristianini N,Bartlett P L,et al.Learning the kernel matrix with semidefinite programming[J].Journal of Machine Learning Research,2004,5:27-72.

[28]Huang H C,Chuang Y Y,Chen C S.Affinity aggregation for spectral clustering[C]//Proceedings of the 2012 IEEE Conference on Computer Vision and Pattern Recognition,Providence,Jun 16-21,2012.Washington:IEEE Computer Society,2012:773-780.

[29]Luxburg U V.A tutorial on spectral clustering[J].Statistics and Computing,2007,17(4):395-416.

[30]Strehl A,Ghosh J.Cluster ensembles-a knowledge reuse framework for combining multiple partitions[J].Journal of Machine Learning Research,2002,3:583-617.