基于改進K-Means聚類醫(yī)學圖像配準

陳 園，侯 贊，劉軍華，雷超陽

（1. 湖南郵電職業(yè)技術(shù)學院互聯(lián)網(wǎng)工程系，湖南 長沙 410015；2. 國網(wǎng)湖南省電力公司檢修公司，湖南 長沙 410004）

0 引言

近年來，針對各種醫(yī)學圖像配準技術(shù)的研究非常熱門，該技術(shù)取得了快速的發(fā)展，醫(yī)療設(shè)備MRI、CT、PET和SPECT等的圖像分辨率有很大的提高，成像質(zhì)量也取得了較好的效果，能為臨床醫(yī)學診斷和治療提供重要可靠的參考依據(jù)。鑒于該技術(shù)的重要性，專家及科研工作者們提出了很多有效解決問題的方法，其中基于圖像特征的配準方法[1-3]得到較為的廣泛應(yīng)用，其原理是通過查找圖像間共有的并有明顯特征的，從而獲取變換參數(shù)。采用該配準方法簡單且大大提高了計算效率。但由于醫(yī)學圖像的復(fù)雜性，圖像特征點的正確提取[4-7]非常重要，采用該方法時計算機很難自動準確提取圖像特征點，需要我們?nèi)斯ぽo助其選取圖像特征點，因此該配準方法在自適應(yīng)性方面還需要進一步提高。

為了解決以上問題，論文在仔細研究 K-Means聚類算法和ICP算法之后，提出了基于改進K-Means聚類醫(yī)學圖像配準（Medical Image Registration Using Improved K-Means Clustering，RIKMC）。該方法通過計算參考圖像和浮動圖像的質(zhì)心，獲得配準平移初始值；對醫(yī)學圖像坐標進行中心化處理，得到兩行坐標矩陣，構(gòu)成KMC（K-Means Clustering）的樣本集合，獲取初始聚類中心，把樣本集合聚成2類；把這 2個聚類中心擬合成一條直線，可以算出該直線的斜率，從而得出其傾斜角，獲得配準旋轉(zhuǎn)初始值；通過其自動選擇特征點，得到參考點集和浮動點集。該方法既可用于單模態(tài)圖像配準，也可以用于多模態(tài)配準。

1 基于改進K-Means聚類醫(yī)學圖像配準

1.1 ICP算法

ICP算法（Iterative Closest Point）最初由Besl和Mckay等人[8-10]提出的（ K =3），假設(shè)在 RK空間存在兩個點集：參考點集 X = { xi,i = 1 ,2,… ,N}且xi=［xi1… xiK］T，和浮動點集 Y = { yj, j = 1 ,2,… ,M }且=…T； Z = {,i=1,2,…,M}為浮動 點集Y在參考點集X得到的最近點集， zi=…］T且∈X；其中T代表K×1的平移矩0陣，R0代表K×K的旋轉(zhuǎn)矩陣。得到最小化目標函數(shù)：

其中，用單位四元組可以計算得到旋轉(zhuǎn)矩陣0R[9-10]，然后計算平移矩陣0T：

從以一描述中可以看到，ICP計算量非常大，受浮動點集初始旋轉(zhuǎn)矩陣和平移矩陣的影響較大，很容易使目標函數(shù)陷入局部極小。

1.2 獲取圖像的質(zhì)心坐標

對于二維離散函數(shù) f (x,y)，它的( p + q)階矩定義為[11-12]：

參數(shù)( )+p q稱為矩的階。

零階矩是物體的面積，其定義為[11-12]：

所有的一階矩和高階矩除以0,0M 后，與物體的大小無關(guān)。

當 1=p ， 0=q 和 0=p ， 1=q 時[11-12]，

則稱(x,y)為圖像中的質(zhì)心坐標。由（3）（4）（5）式計算出參考圖像R、浮動圖像F的零階矩與一階矩，從而得到它們的質(zhì)心坐標為(xR, yR)和(xF, yF)。

假設(shè)(x,y)為物體的質(zhì)心，則中心矩的定義為[11-12]：

1.3 使用改進的KMC獲取圖像的傾斜角

MacQueen J.B在1967年提出了K-means聚類算法[13-15]（K-Means Clustering,KMC），其前提條件是在最小化誤差函數(shù)條件下，將數(shù)據(jù)劃分為給定的聚類數(shù)并能處理大量數(shù)據(jù)。但該算法也存在一些缺陷和不足，那就是在運行前要先指定收斂條件、初始聚類中心、聚類數(shù)和迭代次數(shù)，再按照相似性原則將每個樣本分布到最相似或最近的聚類中心而構(gòu)成類后進行重新分配，通過反復(fù)迭代從而達到類收斂。

1.3.1 KMC算法

設(shè)X ={ x1, x2，…，xn}是一個由n個樣本組成的集合，且有 x =[x ,x ,… ,x ]T(i = 1,2,… ,n)，c是i i,1i,2i,l給定的聚類數(shù)， vj( j = 1,2， … ,c)是每個聚類的中心，且有 vj= [ vj,1, vj,2, … ,vj,l]T，每個聚類集合由 rj個樣本組成， xkj表示第 j類的第k個樣本， j = 1 ,2,… ,c，k =1,2,… ,rj。則每個聚類中心由下式定義：

因此，聚類目標函數(shù)可定義為：

KMC聚類時，一般選擇 d (xi, vj)為歐氏距離，即d()= x -v2。該算法過度依賴于初始聚類j i j中心的選擇，如果選擇不當，使分類的結(jié)果會嚴重偏離全局最優(yōu)分類，特別是聚類數(shù)較大時其偏離最優(yōu)值的缺點更為明顯，就需要通過更多次數(shù)的聚類，得到的效果才有可能較為滿意。

1.3.2 KMC初始聚類中心選擇

由于 KMC對初始聚類中心過度依賴及敏感從而造成聚類結(jié)果不穩(wěn)定。如圖1（a）和圖1（b）兩幅傾斜醫(yī)學圖像，當我們隨機選擇初始聚類中心時，采用KMC方法各自運行10次后，得到的結(jié)果如圖2（a）和圖2（b）所示。

圖1 傾斜醫(yī)學圖像Fig.1 Tilt medical images

從圖2（a）和圖2（b）中可以看到，采用KMC方法獲得的傾斜角在一個范圍內(nèi)震蕩，很不穩(wěn)定并會形成多個極值點，從而得到的傾斜角不可靠，使得聚類結(jié)果不可預(yù)測，這樣醫(yī)學圖像傾斜角的不確定性將嚴重影響其圖像后續(xù)的進一步處理，還有KMC方法運行時間相對波動較大。由于該算法過度依賴于初始聚類中心的選擇，所以我們一定要適當選擇初始聚類中心。

圖2 隨機選取初始聚類中心與得到相應(yīng)的傾斜角、KMC運行時間的關(guān)系Fig.2 The relationship among the random selectionof initial clustering center, the derived angle,and the running time of KMC

在論文中，生成初始聚類中心的方法如下：一是按照樣本集X的大小n計算出 h alf =integer[n / 2]；二是將樣本集分為兩類聚類集合，一聚類集合為前half樣本，另一聚類集合為后(n - half)樣本；三是計算每類的初始聚類中心，如下式：

對于圖 1（a）、圖 1（b）兩幅傾斜醫(yī)學圖像，通過（9）式計算出初始聚類中心，得到如圖3（a）和圖3（b）所示的關(guān)系，從此可以看出，傾斜角較為穩(wěn)定，運行時間波動也不大，這有利于獲得較好的醫(yī)學圖像傾斜角。

圖3 通過（9）式生成初始聚類中心與獲得的傾斜角、KMC運行時間的關(guān)系Fig.3 The relationship among the clustering center produced by Equation (9), the derived angle, and the running time of KMC

1.3.3 獲取圖像的傾斜角

綜上所述，改進 K-Means聚類獲取醫(yī)學圖像傾斜角（Acquisition of Rotation Angle Based on Improved K-Means Clustering，AIKMC）過程描述如下：

步驟1：對于傾斜圖像F的邊緣檢測采用B樣條梯度算法得到二值化的圖像B；

步驟2：找出圖像B包圍盒的四周邊界；

步驟 3：在包圍盒的四周邊界范圍內(nèi)截取相應(yīng)子圖像 F _ Sub；

步驟 4：通過（3）式和（4）式計算出子圖像F _ Sub 的一階矩 M1,0、 M0,1與零階矩 M0,0；

步驟 5：通過（5）式計算出子圖像 F _ Sub的質(zhì)心坐標(x,y)；

步驟6：對子圖像 F _ Sub建立傾斜圖像坐標矩陣PR，其中H×W代表元素個數(shù)（行數(shù)），它們之間的關(guān)系為∈R(H×W)×2：

其中， i = 1,2,… ,H; j = 1,2,… ,W 。

步驟7：根據(jù)傾斜圖像矩陣PR的大小n=H×W，計算 h alf =integer[n / 2]；然后把 PR劃分為兩類：以PR前half個樣本為一聚類集合，后(n - half)樣本為另一聚類集合，根據(jù)（9）式生成KMC初始聚類中心；

步驟8：通過KMC，把 PR聚成2類 vj( j =1,2)，且]T；

步驟 9：把 2個聚類中心擬合成一條直線，并計算該直線斜率k；

步驟10：根據(jù)k，得到傾斜角α：

1.4 基于改進K-Means聚類醫(yī)學圖像配準

由于傳統(tǒng)的ICP算法存在上述問題，論文提出了基于改進 K-Means聚類醫(yī)學圖像配準（Medical Image Registration Using Improved K-Means Clustering，RIKMC）方法，具體描述如下：

步驟1：根據(jù)圖像矩和AIKMC，獲取參考圖像R和浮動圖像F的質(zhì)心坐標(xr, yr)、(xf, yf)，以及旋轉(zhuǎn)角αr和αf；

步驟 2：計算配準初始值：Δx = ( xf- yf)，Δ y = yf-yr， Δα = αf- αr；

步驟3：把Δx、Δy和Δα作為ICP算法的初始值，即：

步驟 4：對得到的參考圖像R與浮動圖像F進行邊緣檢測得到二值化圖像RB和FB，定義灰度值為 0和1；

步驟5：將灰度值為1的二值化圖像RB 和FB 的像素點的坐標看做是ICP算法的參考圖像點集X和浮動圖像點集Y；

步驟 6：通過 ICP方法進行迭代后得到所需要的平移矩陣和旋轉(zhuǎn)矩陣。

2 實驗?zāi)M及結(jié)果分析

實驗?zāi)M軟件為Matlab7.1，操作系統(tǒng)Windows XP，內(nèi)存2 GB，CPU為Pentium? Dual-Core E5500 2.80 GHz，實驗數(shù)據(jù)選取了美國 Vanderbilt大學Retrospective Registration Evaluation Projection（RREP）項目組的國際通用剛性配準腦部圖像數(shù)據(jù)。誤差ρ定義為：

（13）式中，通過配準獲取的浮動圖像變換參數(shù)用Δi表示，浮動圖像相對于參考圖像的標準變換參數(shù)用Δsi表示。得到總誤差ρ，即：

2.1 單模態(tài)醫(yī)學圖像配準

在該實驗中參考圖像選取PET第二層圖像、頭部CT patient_007第三層圖像和MR_PD_rectified第四層圖像，它們的圖像灰度都為256級，其中PET圖像分辨率像素為128128×，CT圖像分辨率像素為512512×，MR圖像分辨率像素為256256×，如參考圖像圖4（a）、5（a）、6（a）所示。通過表1單模態(tài)配準時浮動圖像實際變換參數(shù)，將相對應(yīng)的參考圖像通過平移、旋轉(zhuǎn)后得到相對應(yīng)的浮動圖像，如浮動圖像圖4（b）和4（c）、圖 5（b）和5（c）、圖6（b）和6（c）所示。

我們首先使用 ICP來完成單模態(tài)配準，然后使用RIKMC來實現(xiàn)單模態(tài)配準。實驗圖像如圖4、圖5和圖6所示，實驗結(jié)果如圖7、圖8以及表2所示。

從表2可以得出如下結(jié)果，在單模態(tài)配準時該方法相對于ICP配準方法有著明顯的優(yōu)勢且整個配準時間較短；特別是圖像較大時，RIKMC配準速度優(yōu)勢更明顯，但是當圖像較小，RIKMC和ICP配準速度相差不大。至于配準精度，ICP配準完全失敗，平移變換陷入局部最優(yōu)；而RIKMC均能成功配準，精度較高，這也證明了ICP配準成功與否，依賴于初始值的選擇，而論文方法使用圖像慣量矩陣產(chǎn)生的初始值，已經(jīng)接近最優(yōu)值，不但節(jié)省了尋優(yōu)時間，而且避免陷入局部最優(yōu)。具體說來，RIKMC對CT1、CT2、MR1、MR2圖像配準成功，精度較高，但對PET1、PET2雖然能成功配準，但是配準精度有待提高，特別是平移參數(shù)誤差較大。

圖4 圖像CTFig.4 CT experimental images

圖5 圖像MRFig.5 MR experimental images

圖6 圖像PETFig.6 PET experimental images

表1 單模態(tài)醫(yī)學圖像配準時浮動圖像實際變換參數(shù)Tab.1 Transformation parameters of mono-modality floating images

圖7 ICPFig.7 ICP

圖8 RIKMCFig.8 RIKMC

表2 單模態(tài)醫(yī)學圖像配準浮動圖像性能比較Tab.2 Performance of registering the mono-modality images

2.2 多模態(tài)醫(yī)學圖像配準

在該實驗中參考圖像選取PET第三層圖像、頭部CT patient_007第五層圖像、MR_PD_rectified第六層圖像。實驗圖像由四組組成：1組如圖9所示，選取參考圖像CT1，浮動圖像MR1，圖像分辨率像素為256256×；2組如圖10所示，選取參考圖像MR2，浮動圖像CT2，圖像分辨率像素為256256×；3組如圖 11所示，選取參考圖像 CT3，浮動圖像PET1，圖像分辨率像素為128128×；4組如圖12所示，選取參考圖像 MR3，浮動圖像 PET2，圖像分辨率像素為128128×，所有圖像灰度均為256級。每幅參考圖像變換后的浮動圖像的實際變換參數(shù)如表3所示。

圖9 1組Fig.9 The first group of images

圖10 2組Fig.10 The second group of images

圖11 3組Fig.11 The third group of images

圖12 4組Fig.12 The fourth group of images

表3 多模態(tài)醫(yī)學圖像配準時浮動圖像實際變換參數(shù)比較Tab.3 Transformation parameters of multi-modality floating images

我們首先使用ICP來完成單模態(tài)配準，然后使用RIKMC來實現(xiàn)單模態(tài)配準。實驗結(jié)果如圖13、圖14、圖15以及表4所示。

從表4可以看到，在多模態(tài)配準時，基于ICP的配準方法和RIKMC配準方法相比較，RIKMC配準方法整個配準時間較短；特別是圖像較大時，RIKMC配準速度優(yōu)勢更明顯，但是當圖像較小，RIKMC和ICP配準速度相差不大。至于配準精度，ICP對2組MR2和CT2配準基本成功，對于1組CT1和MR1、3組 CT3和PET1、4組MR3和PET2配準配準完全失敗，陷入局部最優(yōu)，因此整體上講，ICP配準失敗率較高；RIKMC均能成功配準，雖然對3組CT3和PET1能成功配準，但是配準誤差相對較高。因此，從整體上講，MICP也適合多模態(tài)配準。

圖13 1組配準結(jié)果Fig.13 Result figures of the first group

圖14 2組配準結(jié)果Fig.14 Result figures of the second group

圖15 3組配準結(jié)果Fig.15 Result figures of the third group

圖16 4組配準結(jié)果Fig.16 Result figures of the fourth group

表4 多模態(tài)醫(yī)學圖像配準浮動圖像性能比較Tab.4 Performance of registering the multi-modality images

3 結(jié)論

為了解決ICP算法配準速度慢、配準成功率低等問題，論文利用改進的K-Means聚類方法，提出了基于改進K-Means聚類醫(yī)學圖像配準算法。通過實驗得出該算法既可用于單模態(tài)圖像配準，也可以用于多模態(tài)圖像配準；具有運算量少、圖像配準速度較快、計算比較簡單、精確度較高等特點，并且解決了圖像配準容易陷入局部最優(yōu)的問題。

[1] Nejati M, Pourghassem H. Multiresolution Image Registration in Digital X-Ray Angiography with Intensity Variation Modeling[J]. Journal of Medical Systems, 2014, 38(1): 10-18.

[2] Pan MS, Jiang JJ, Rong QS. A modified medical image registration[J]. Multimedia Tools and Applications, 2014,70(3): 1585–1615.

[3] J B Antoine Maintz，M A Viergever．A survey of medical image registration[J]. Medical Image Analysis, 1998, 2(1): 1-36.[4] Hyunjin P, Peyton H, Kristy K, et a1. Adaptive registration using local information measures[J]. Medical Image Analysis,2004, 8(4): 465-473.

[5] Can A, Stewart C. A feature-based, robust, hierarchical algorithm for registration palm of images of the curved human retina[J]．IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 24(3): 347-363.

[6] Liu Y, Cheng HD, Huang JH. An Effective Non-rigid Registration Approach for Ultrasound Image Based On“Demons” Algorithm[J]. Journal of Digital Imaging, 2013,26(3): 521-529.

[7] Liu W, Ribeiro E. Incremental variations of image moments for nonlinear image registration[J]. Signal, Image and Video Processing, 2014, 8(3): 423-432.

[8] Arun KS, Huang TS, Blostein SD. Least-squares fitting of two 3-D point sets[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1987, 9(5): 698-700.

[9] Besl PJ, Mckay ND. A method for registration of 3-D shapes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14(2): 239-256.

[10] Kaneko S., Kondo T., Miyamoto A.Robust matching of 3D contours using iterative closest point algorithm improved by M-estimation[J]. Pattern Recognition, 2003, 36(9): 2041-2047.

[11] Hu MK. Visual pattern recognition by moment invariants[J].IRE Transactions On Information Theory, 1962, 8(2): 179-187.

[12] Wong YR. Matching with invariant moments[J]. Computer Graphics and Image Processing, 1978, 8(1): 16-24.

[13] Anil K J. Data clustering: 50 years beyond K-Means[J].Pattern Recognition Letters, 2010, 31(8): 651-666.

[14] Mahajan M, Nimbor P, Varadarajan K. The planar K-means problem is NP-hard[J]. Lecture Notes in Computer Science,2009(5431): 274-285.

[15] Lai J Z C, Tsung-Jen H. Fast global K -means clustering using cluster membership and inequality[J]. Pattern Recognition,2010(43): 1954-1963.