王保加，潘海為，謝曉芹，張志強(qiáng)，馮曉寧

哈爾濱工程大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，哈爾濱 150001

1 引言

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

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

然而醫(yī)學(xué)圖像特征提取和普通圖像相比有很大的難點(diǎn)，醫(yī)學(xué)圖像的相似性度量涉及到病理學(xué)和CT成像相關(guān)知識(shí)，導(dǎo)致普通圖像的特征提取方法無(wú)法完全表達(dá)醫(yī)學(xué)圖像的全部信息。比如對(duì)于在不同區(qū)域發(fā)生相同病變的醫(yī)學(xué)圖像，它們?cè)诓±韺W(xué)上應(yīng)歸屬同類別的病變圖像，但現(xiàn)有的特征提取方法忽略了這種高層語(yǔ)義信息，不僅不能突顯不同位置尤其是病變區(qū)域的重要性，還會(huì)因?yàn)檎^(qū)域和病變區(qū)域的相互干擾造成圖像相似性比較不準(zhǔn)確。同時(shí)醫(yī)生對(duì)CT圖像診斷時(shí)，也更加關(guān)注具有診斷意義的ROI（region of interest）區(qū)域，而忽略細(xì)微紋理的差異。因此提取出能夠反映人類視覺(jué)認(rèn)知的醫(yī)學(xué)ROI特征，在醫(yī)學(xué)圖像的聚類中具有重要的意義。

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

針對(duì)以上研究存在的問(wèn)題，本文提出了一種融合圖像紋理特征和醫(yī)學(xué)圖像特有領(lǐng)域知識(shí)的形態(tài)學(xué)特征的多模態(tài)特征抽取方法。首先利用紋理特征在相似紋理圖像上的良好區(qū)分效果，進(jìn)而融入醫(yī)學(xué)圖像領(lǐng)域知識(shí)，結(jié)合醫(yī)學(xué)圖像ROI在病變圖像的區(qū)分能力，實(shí)現(xiàn)了對(duì)復(fù)雜腦部CT圖像集的聚類，并使用了多核譜聚類算法驗(yàn)證了聚類的效果。

本文組織結(jié)構(gòu)如下：第2章介紹多模態(tài)特征醫(yī)學(xué)圖像聚類的總體框架；第3章研究紋理特征和形態(tài)學(xué)特征提取方法；第4章討論基于多模態(tài)特征的聚類方法；第5章展示實(shí)驗(yàn)結(jié)果和相關(guān)分析；最后總結(jié)全文。

2 多模態(tài)特征聚類總體框架

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

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

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

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

3 特征提取

3.1 紋理特征

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

3.1.1 紋理特征描述

（1）LBP

LBP特征通過(guò)比較中心像素點(diǎn)和鄰域半徑內(nèi)像素點(diǎn)的灰度值，計(jì)算出中心像素點(diǎn)二進(jìn)制編碼序列，再統(tǒng)計(jì)不同LBP序列的數(shù)目來(lái)表達(dá)圖像的LBP紋理特征FL。對(duì)于一個(gè)像素點(diǎn)的LBP值定義為：

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

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

Fig.3 Calculate LBP圖3 LBP值計(jì)算

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

（2）GLCM

GLCM是像素距離和角度的矩陣。本文使用灰度共生矩陣的角二階矩、對(duì)比度、相關(guān)性和熵值來(lái)提取GLCM的紋理特征

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

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

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

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

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

例2分別從 0°、45°、90°、135°共4個(gè)方向計(jì)算出圖4（a）的GLCM的16個(gè)相關(guān)值，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的提取過(guò)程。首先將圖像4（a）進(jìn)行canny邊緣檢測(cè)，得到邊緣紋理圖像4（b）；然后對(duì)邊緣紋理上的像素點(diǎn)I分別沿x軸、y軸計(jì)算邊緣方向梯度值dx和dy，則像素點(diǎn)I的邊緣方向如圖 4（c）對(duì) [-180°,180°]角度區(qū)間每隔 10°劃分一組，一共分成36組邊緣方向區(qū)間，統(tǒng)計(jì)像素點(diǎn)落在每組區(qū)間的個(gè)數(shù)，最終生成36維邊緣方向向量FE。

Fig.4 EDH extraction圖4EDH提取圖

3.1.2 FET方法

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

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

算法1FET算法

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

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

1.ForIM的每一張圖像IMi

2.ForIMi的每一個(gè)分割區(qū)域j

3.計(jì)算第j區(qū)域的添加到向量數(shù)組

4.End for

6.End for

7.ForXD×Z的每一列Xi

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

9.End for

10.計(jì)算協(xié)方差矩陣：

11.求出協(xié)方差矩陣Z個(gè)最大的特征值及對(duì)應(yīng)的特征向量UZ×D=[μ1,μ2,…,μZ]T；

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

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

3.2 形態(tài)學(xué)特征

醫(yī)學(xué)圖像ROI的提取是為了獲取醫(yī)學(xué)圖像中具有診斷意義的像素區(qū)域。本節(jié)給出形態(tài)特征ROI提取的方法、ROI相關(guān)描述和圖像ROI向量化的方法。

3.2.1 REGA方法

本文使用具有醫(yī)學(xué)領(lǐng)域知識(shí)的黑影和白影區(qū)域作為ROI區(qū)域。以下是ROI提取的過(guò)程，定義1給出了一個(gè)圖像集中ROI的相關(guān)概念。

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

本文提出了ROI區(qū)域提取算法REGA，算法的基本思想是根據(jù)黑影灰度閥值BGrayε、白影灰度閥值WGrayε和面積閥值A(chǔ)reaε的限制分割出最佳的ROI區(qū)域，具體過(guò)程如算法2所示。

算法2REGA算法

輸入：圖像集IM中的第i張圖像IMi；灰度閥值BGrayε和WGrayε；面積閥值A(chǔ)reaε。

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

1.IMi使用sobel算子進(jìn)行邊緣檢測(cè)、高斯濾波噪聲抑制、形態(tài)學(xué)開(kāi)閉操作平滑去噪；

2.IMi提取區(qū)域邊緣并進(jìn)行二值化處理，得到分割圖像IMB；

3.For圖像IMB分割區(qū)域

4.IfIMB分割區(qū)域面積＜Areaε

5. ThenIMB分割區(qū)域值設(shè)置為0；

6.End if

7.Else ifIMi邊緣灰度值＜BGrayε或邊緣灰度值＞W(wǎng)Grayε

8. ThenIMB分割區(qū)域值設(shè)置為1；

9. ElseIMB分割區(qū)域值設(shè)置為0；

10. End if

11.End for；

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

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

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

對(duì)于提取到的ROI，從形態(tài)學(xué)的特征上可以很好地描述不同ROI的特征。本文使用ROI的面積、ROI的位置、ROI的離心率、ROI的周徑比、ROI的圓形度、ROI的方向和ROI黑影白影標(biāo)記來(lái)表達(dá)一幅圖像。對(duì)于一個(gè)ROI區(qū)域R，具體定義如下：

R的面積（Area）表示為AreaR=M，其中M為ROI區(qū)域的像素點(diǎn)的總個(gè)數(shù)。

R的重心位置（Location）表示為L(zhǎng)ocR=(Centroidx,Centroidy)，其中Centroidx、Centroidy表示該區(qū)域的二維重心坐標(biāo)值。

Fig.5 ROI area extraction圖5 ROI區(qū)域提取圖

R的離心率（Eccentricity）表示為其中LAL表示與R具有相同標(biāo)準(zhǔn)二階中心矩的橢圓的長(zhǎng)軸長(zhǎng)度（long axis length，LAL），SAL表示與R具有相同標(biāo)準(zhǔn)二階中心矩的橢圓的短軸長(zhǎng)度（short axis length，SAL）。

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

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

R的方向（Orientation）表示為OriR=θ，其中θ為與R具有相同標(biāo)準(zhǔn)二階中心矩的橢圓的長(zhǎng)軸與X軸正方向的夾角。

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

因此任意的兩個(gè)ROI區(qū)域Ri可由一個(gè)k=7維向量和灰度標(biāo)記表示。對(duì)于任意的兩個(gè)ROI區(qū)域Ri和Rj，距離計(jì)算為：

3.2.2 ROIV方法

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

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

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

定義3R灰度級(jí)（gray level，GL）為R內(nèi)所有像素點(diǎn)灰度的平均值，其中M為R中像素點(diǎn)的總數(shù)，Ii表示第i個(gè)像素點(diǎn)的灰度值。

定義4圖像集IM的正?；叶燃?jí)（normal gray Level，NGL）為正常圖像樣本集ROI灰度級(jí)的平均值，，其中SR為正常圖像提取到的ROI數(shù)目，GLi為第i個(gè)ROI區(qū)域Ri的灰度級(jí)。

定義5對(duì)于一個(gè)ROI區(qū)域R，R的灰度級(jí)權(quán)重（gray level weight，GLW）為：

定義6對(duì)于圖像集IM中圖像IMi的第j個(gè)ROI區(qū)域Rij，Rij面積權(quán)重（area weight，AW）為：

RTW-IDF加權(quán)：

RTW：對(duì)于圖像集IM中圖像IMi的第j個(gè)ROI區(qū)域Rij，Rij權(quán)重（ROI weight，RW）為：

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

圖像IMi的向量Vi的第j維中包含的單詞權(quán)值（term weight，TW）和為為聚類中IMi在第j類中RT數(shù)目，RWj為第j個(gè)RT權(quán)值。

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

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

算法3ROIV算法

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

輸出：圖像集IM的形態(tài)學(xué)特征矩陣。

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

2.利用式（7）對(duì)RD進(jìn)行K-Means聚成C類；

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

4.ForIMi的每一張圖像

5.ForIMi的每一個(gè)ROI區(qū)域Rij

6. 歐氏距離計(jì)算Rij最近U的聚類中心，標(biāo)記為類別p添加數(shù)組q，式（9）計(jì)算Rij在圖像IMi中的ROI權(quán)值rw；

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

8.End for

9.For數(shù)組q的每一個(gè)元素q[i]

10.取回類別標(biāo)記p=q[i]；

計(jì)算MatrixRip+=w[i]；

11.End for

12.End for

13.返回MatrixRN×C；

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

4 基于多模態(tài)特征的聚類算法

4.1 MSF方法

相似性度量選取對(duì)聚類結(jié)果有很大的影響，高斯核函數(shù)常被用于度量?jī)蓮垐D像的相似度：

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

其中，σi表示數(shù)據(jù)xi與第L(L=7)個(gè)近鄰點(diǎn)的距離。

第3章利用特征提取算法得到了圖像集IM的紋理特征向量矩陣MatrixT和形態(tài)學(xué)特征向量矩陣MatrixB?？紤]到兩種特征在圖像相似度上的互補(bǔ)作用，本文重新定義了多模態(tài)相似性度量函數(shù)MSF：

其中，φi表示圖像IMi分別使用紋理特征和形態(tài)學(xué)特征前L個(gè)近鄰點(diǎn)集合中公共交集的圖像數(shù)目。設(shè)A表示IMi在紋理特征度量下前L個(gè)近鄰圖像的集合，B表示IMi在形態(tài)學(xué)特征度量下前L個(gè)近鄰圖像的集合，則兩個(gè)集合中相同的圖片數(shù)量φi=A?B。

通過(guò)本文定義的MSF計(jì)算紋理特征向量矩陣MatrixT得到N×N的紋理特征相似度對(duì)稱矩陣同理計(jì)算形態(tài)學(xué)特征向量矩陣MatrixR得到形態(tài)學(xué)特征相似度對(duì)稱矩陣SMatrixR。

4.2 MKSC方法

多核學(xué)習(xí)（multi-kernel learning，MKL）是由Lanckriet等人提出的[27]，在此基礎(chǔ)上Huang等人又把多核學(xué)習(xí)方法擴(kuò)展到非監(jiān)督的聚類領(lǐng)域[28]，多核學(xué)習(xí)可以自適應(yīng)地完成對(duì)不同特征分配權(quán)值。

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

其中，Wij表示數(shù)據(jù)xi和xj的相似度。

對(duì)于4.1節(jié)中數(shù)據(jù)集兩種特征相似度矩陣，給定紋理特征權(quán)值為α1，形態(tài)學(xué)特征權(quán)值為α2，且α1+α2=1，則式（12）轉(zhuǎn)變?yōu)榍蠼鈳Ъs束條件的式（13）：

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

算法4MKSC算法

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

輸出：K個(gè)聚類類別 {Y1,Y2,…YK}，特征權(quán)值α1、α2。

1.初始化特征權(quán)值α1=α2=0.5；

2.Repeat

3.構(gòu)造相似性度量矩陣：

5.計(jì)算L最小的K個(gè)特征值λ1,λ2,…,λK所對(duì)應(yīng)的特征列向量u1,u2,…,uK；

7.更新α1、α2；

8.Until收斂

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

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

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

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

5 實(shí)驗(yàn)

5.1 數(shù)據(jù)集和環(huán)境

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

5.2 評(píng)價(jià)標(biāo)準(zhǔn)

本文采用標(biāo)準(zhǔn)化互信息（normalized mutual information，NMI）作為評(píng)價(jià)標(biāo)準(zhǔn)[30]。對(duì)于具有K個(gè)真實(shí)分類集T和聚類結(jié)果集C的NMI：

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

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

5.3 實(shí)驗(yàn)設(shè)置和結(jié)果

紋理特征提取時(shí)需要對(duì)圖像進(jìn)行分割，不同的分割方式會(huì)有不同的結(jié)果，通常分割后圖像的總數(shù)小于10[25]?；诖吮緦?shí)驗(yàn)比較了5種不同的圖像分塊分割方式：1×1（不分區(qū)），2×2，3×3，2×4和4×2。表1為不同分塊模式下的NMI值，當(dāng)分割為4×2的模式時(shí)，能取得最好的效果，因此本實(shí)驗(yàn)使用4×2的模式進(jìn)行圖像分塊。

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

本文還需要對(duì)圖像ROI向量化時(shí)向量維數(shù)C值進(jìn)行設(shè)置，太小的C值會(huì)把不同類別的ROI聚成一類，太大的C值會(huì)使聚類太分散，造成圖像區(qū)分不明顯。圖6為不同C值下使用形態(tài)學(xué)特征圖像聚類的NMI值。當(dāng)C取值在28～31附近時(shí)NMI值趨于最大且穩(wěn)定于0.58，因此本文以下實(shí)驗(yàn)C都設(shè)置為30。

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

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

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

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

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

本文提出的醫(yī)學(xué)圖像聚類方法不需要復(fù)雜的訓(xùn)練過(guò)程，實(shí)驗(yàn)的NMI值雖然在低數(shù)量級(jí)下低于Zhan（2016）的方法，但隨著圖像數(shù)目的增加實(shí)驗(yàn)結(jié)果呈現(xiàn)增長(zhǎng)趨勢(shì)。這是由于圖像集越大，醫(yī)學(xué)圖像的ROI類型越多越全面，對(duì)醫(yī)學(xué)圖像的特征描述更準(zhǔn)確。本實(shí)驗(yàn)的方法整體優(yōu)于Liang（2016）方法。

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

6 結(jié)束語(yǔ)

本文提出了一種融合圖像紋理特征和醫(yī)學(xué)圖像特有領(lǐng)域知識(shí)的形態(tài)學(xué)特征的多模態(tài)特征抽取方法，旨在解決傳統(tǒng)特征提取方法不能有效地表達(dá)醫(yī)學(xué)圖像特征導(dǎo)致聚類不準(zhǔn)確的難題。本文提出了醫(yī)學(xué)圖像的ROI提取算法和ROI形態(tài)學(xué)特征描述，使用特征融合的多核譜聚類算法完成圖像聚類。實(shí)驗(yàn)結(jié)果表明此方法可以減小非正常醫(yī)學(xué)圖像和正常醫(yī)學(xué)圖像聚類的相互干擾，實(shí)現(xiàn)復(fù)雜醫(yī)學(xué)圖像聚類。

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

[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.