陳 園，侯 贊，劉軍華，雷超陽(yáng)

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

0 引言

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

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

1 基于改進(jìn)K-Means聚類醫(yī)學(xué)圖像配準(zhǔn)

1.1 ICP算法

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

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

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

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

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

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

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

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

當(dāng) 1=p ， 0=q 和 0=p ， 1=q 時(shí)[11-12]，

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

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

1.3 使用改進(jìn)的KMC獲取圖像的傾斜角

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

1.3.1 KMC算法

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

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

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

1.3.2 KMC初始聚類中心選擇

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

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

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

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

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

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

圖3 通過（9）式生成初始聚類中心與獲得的傾斜角、KMC運(yùn)行時(shí)間的關(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 獲取圖像的傾斜角

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

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

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

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

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

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

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

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

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

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

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

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

1.4 基于改進(jìn)K-Means聚類醫(yī)學(xué)圖像配準(zhǔn)

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

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

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

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

步驟 4：對(duì)得到的參考圖像R與浮動(dòng)圖像F進(jìn)行邊緣檢測(cè)得到二值化圖像RB和FB，定義灰度值為 0和1；

步驟5：將灰度值為1的二值化圖像RB 和FB 的像素點(diǎn)的坐標(biāo)看做是ICP算法的參考圖像點(diǎn)集X和浮動(dòng)圖像點(diǎn)集Y；

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

2 實(shí)驗(yàn)?zāi)M及結(jié)果分析

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

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

2.1 單模態(tài)醫(yī)學(xué)圖像配準(zhǔn)

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

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

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

圖4 圖像CTFig.4 CT experimental images

圖5 圖像MRFig.5 MR experimental images

圖6 圖像PETFig.6 PET experimental images

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

圖7 ICPFig.7 ICP

圖8 RIKMCFig.8 RIKMC

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

2.2 多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)

在該實(shí)驗(yàn)中參考圖像選取PET第三層圖像、頭部CT patient_007第五層圖像、MR_PD_rectified第六層圖像。實(shí)驗(yàn)圖像由四組組成：1組如圖9所示，選取參考圖像CT1，浮動(dòng)圖像MR1，圖像分辨率像素為256256×；2組如圖10所示，選取參考圖像MR2，浮動(dòng)圖像CT2，圖像分辨率像素為256256×；3組如圖 11所示，選取參考圖像 CT3，浮動(dòng)圖像PET1，圖像分辨率像素為128128×；4組如圖12所示，選取參考圖像 MR3，浮動(dòng)圖像 PET2，圖像分辨率像素為128128×，所有圖像灰度均為256級(jí)。每幅參考圖像變換后的浮動(dòng)圖像的實(shí)際變換參數(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ī)學(xué)圖像配準(zhǔn)時(shí)浮動(dòng)圖像實(shí)際變換參數(shù)比較Tab.3 Transformation parameters of multi-modality floating images

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

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

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

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

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

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

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

3 結(jié)論

為了解決ICP算法配準(zhǔn)速度慢、配準(zhǔn)成功率低等問題，論文利用改進(jìn)的K-Means聚類方法，提出了基于改進(jìn)K-Means聚類醫(yī)學(xué)圖像配準(zhǔn)算法。通過實(shí)驗(yàn)得出該算法既可用于單模態(tài)圖像配準(zhǔn)，也可以用于多模態(tài)圖像配準(zhǔn)；具有運(yùn)算量少、圖像配準(zhǔn)速度較快、計(jì)算比較簡(jiǎn)單、精確度較高等特點(diǎn)，并且解決了圖像配準(zhǔn)容易陷入局部最優(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.