基于改进K-Means聚类医学图像配准

陈 园，侯 赞，刘军华，雷超阳

（1. 湖南邮电职业技术学院互联网工程系，湖南 长沙 410015；2. 国网湖南省电力公司检修公司，湖南 长沙 410004）

0 引言

近年来，针对各种医学图像配准技术的研究非常热门，该技术取得了快速的发展，医疗设备MRI、CT、PET和SPECT等的图像分辨率有很大的提高，成像质量也取得了较好的效果，能为临床医学诊断和治疗提供重要可靠的参考依据。鉴于该技术的重要性，专家及科研工作者们提出了很多有效解决问题的方法，其中基于图像特征的配准方法[1-3]得到较为的广泛应用，其原理是通过查找图像间共有的并有明显特征的，从而获取变换参数。采用该配准方法简单且大大提高了计算效率。但由于医学图像的复杂性，图像特征点的正确提取[4-7]非常重要，采用该方法时计算机很难自动准确提取图像特征点，需要我们人工辅助其选取图像特征点，因此该配准方法在自适应性方面还需要进一步提高。

为了解决以上问题，论文在仔细研究 K-Means聚类算法和ICP算法之后，提出了基于改进K-Means聚类医学图像配准（Medical Image Registration Using Improved K-Means Clustering，RIKMC）。该方法通过计算参考图像和浮动图像的质心，获得配准平移初始值；对医学图像坐标进行中心化处理，得到两行坐标矩阵，构成KMC（K-Means Clustering）的样本集合，获取初始聚类中心，把样本集合聚成2类；把这 2个聚类中心拟合成一条直线，可以算出该直线的斜率，从而得出其倾斜角，获得配准旋转初始值；通过其自动选择特征点，得到参考点集和浮动点集。该方法既可用于单模态图像配准，也可以用于多模态配准。

1 基于改进K-Means聚类医学图像配准

1.1 ICP算法

ICP算法（Iterative Closest Point）最初由Besl和Mckay等人[8-10]提出的（ K =3），假设在 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的旋转矩阵。得到最小化目标函数：

其中，用单位四元组可以计算得到旋转矩阵0R[9-10]，然后计算平移矩阵0T：

从以一描述中可以看到，ICP计算量非常大，受浮动点集初始旋转矩阵和平移矩阵的影响较大，很容易使目标函数陷入局部极小。

1.2 获取图像的质心坐标

对于二维离散函数 f (x,y)，它的( p + q)阶矩定义为[11-12]：

参数( )+p q称为矩的阶。

零阶矩是物体的面积，其定义为[11-12]：

所有的一阶矩和高阶矩除以0,0M 后，与物体的大小无关。

当 1=p ， 0=q 和 0=p ， 1=q 时[11-12]，

则称(x,y)为图像中的质心坐标。由（3）（4）（5）式计算出参考图像R、浮动图像F的零阶矩与一阶矩，从而得到它们的质心坐标为(xR, yR)和(xF, yF)。

假设(x,y)为物体的质心，则中心矩的定义为[11-12]：

1.3 使用改进的KMC获取图像的倾斜角

MacQueen J.B在1967年提出了K-means聚类算法[13-15]（K-Means Clustering,KMC），其前提条件是在最小化误差函数条件下，将数据划分为给定的聚类数并能处理大量数据。但该算法也存在一些缺陷和不足，那就是在运行前要先指定收敛条件、初始聚类中心、聚类数和迭代次数，再按照相似性原则将每个样本分布到最相似或最近的聚类中心而构成类后进行重新分配，通过反复迭代从而达到类收敛。

1.3.1 KMC算法

设X ={ x1, x2，…，xn}是一个由n个样本组成的集合，且有 x =[x ,x ,… ,x ]T(i = 1,2,… ,n)，c是i i,1i,2i,l给定的聚类数， 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。则每个聚类中心由下式定义：

因此，聚类目标函数可定义为：

KMC聚类时，一般选择 d (xi, vj)为欧氏距离，即d()= x -v2。该算法过度依赖于初始聚类j i j中心的选择，如果选择不当，使分类的结果会严重偏离全局最优分类，特别是聚类数较大时其偏离最优值的缺点更为明显，就需要通过更多次数的聚类，得到的效果才有可能较为满意。

1.3.2 KMC初始聚类中心选择

由于 KMC对初始聚类中心过度依赖及敏感从而造成聚类结果不稳定。如图1（a）和图1（b）两幅倾斜医学图像，当我们随机选择初始聚类中心时，采用KMC方法各自运行10次后，得到的结果如图2（a）和图2（b）所示。

图1 倾斜医学图像Fig.1 Tilt medical images

从图2（a）和图2（b）中可以看到，采用KMC方法获得的倾斜角在一个范围内震荡，很不稳定并会形成多个极值点，从而得到的倾斜角不可靠，使得聚类结果不可预测，这样医学图像倾斜角的不确定性将严重影响其图像后续的进一步处理，还有KMC方法运行时间相对波动较大。由于该算法过度依赖于初始聚类中心的选择，所以我们一定要适当选择初始聚类中心。

图2 随机选取初始聚类中心与得到相应的倾斜角、KMC运行时间的关系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）两幅倾斜医学图像，通过（9）式计算出初始聚类中心，得到如图3（a）和图3（b）所示的关系，从此可以看出，倾斜角较为稳定，运行时间波动也不大，这有利于获得较好的医学图像倾斜角。

图3 通过（9）式生成初始聚类中心与获得的倾斜角、KMC运行时间的关系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聚类获取医学图像倾斜角（Acquisition of Rotation Angle Based on Improved K-Means Clustering，AIKMC）过程描述如下：

步骤1：对于倾斜图像F的边缘检测采用B样条梯度算法得到二值化的图像B；

步骤2：找出图像B包围盒的四周边界；

步骤 3：在包围盒的四周边界范围内截取相应子图像 F _ Sub；

步骤 4：通过（3）式和（4）式计算出子图像F _ Sub 的一阶矩 M1,0、 M0,1与零阶矩 M0,0；

步骤 5：通过（5）式计算出子图像 F _ Sub的质心坐标(x,y)；

步骤6：对子图像 F _ Sub建立倾斜图像坐标矩阵PR，其中H×W代表元素个数（行数），它们之间的关系为∈R(H×W)×2：

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

步骤7：根据倾斜图像矩阵PR的大小n=H×W，计算 h alf =integer[n / 2]；然后把 PR划分为两类：以PR前half个样本为一聚类集合，后(n - half)样本为另一聚类集合，根据（9）式生成KMC初始聚类中心；

步骤8：通过KMC，把 PR聚成2类 vj( j =1,2)，且]T；

步骤 9：把 2个聚类中心拟合成一条直线，并计算该直线斜率k；

步骤10：根据k，得到倾斜角α：

1.4 基于改进K-Means聚类医学图像配准

由于传统的ICP算法存在上述问题，论文提出了基于改进 K-Means聚类医学图像配准（Medical Image Registration Using Improved K-Means Clustering，RIKMC）方法，具体描述如下：

步骤1：根据图像矩和AIKMC，获取参考图像R和浮动图像F的质心坐标(xr, yr)、(xf, yf)，以及旋转角α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方法进行迭代后得到所需要的平移矩阵和旋转矩阵。

2 实验模拟及结果分析

实验模拟软件为Matlab7.1，操作系统Windows XP，内存2 GB，CPU为Pentium® Dual-Core E5500 2.80 GHz，实验数据选取了美国 Vanderbilt大学Retrospective Registration Evaluation Projection（RREP）项目组的国际通用刚性配准脑部图像数据。误差ρ定义为：

（13）式中，通过配准获取的浮动图像变换参数用Δi表示，浮动图像相对于参考图像的标准变换参数用Δsi表示。得到总误差ρ，即：

2.1 单模态医学图像配准

在该实验中参考图像选取PET第二层图像、头部CT patient_007第三层图像和MR_PD_rectified第四层图像，它们的图像灰度都为256级，其中PET图像分辨率像素为128128×，CT图像分辨率像素为512512×，MR图像分辨率像素为256256×，如参考图像图4（a）、5（a）、6（a）所示。通过表1单模态配准时浮动图像实际变换参数，将相对应的参考图像通过平移、旋转后得到相对应的浮动图像，如浮动图像图4（b）和4（c）、图 5（b）和5（c）、图6（b）和6（c）所示。

我们首先使用 ICP来完成单模态配准，然后使用RIKMC来实现单模态配准。实验图像如图4、图5和图6所示，实验结果如图7、图8以及表2所示。

从表2可以得出如下结果，在单模态配准时该方法相对于ICP配准方法有着明显的优势且整个配准时间较短；特别是图像较大时，RIKMC配准速度优势更明显，但是当图像较小，RIKMC和ICP配准速度相差不大。至于配准精度，ICP配准完全失败，平移变换陷入局部最优；而RIKMC均能成功配准，精度较高，这也证明了ICP配准成功与否，依赖于初始值的选择，而论文方法使用图像惯量矩阵产生的初始值，已经接近最优值，不但节省了寻优时间，而且避免陷入局部最优。具体说来，RIKMC对CT1、CT2、MR1、MR2图像配准成功，精度较高，但对PET1、PET2虽然能成功配准，但是配准精度有待提高，特别是平移参数误差较大。

图4 图像CTFig.4 CT experimental images

图5 图像MRFig.5 MR experimental images

图6 图像PETFig.6 PET experimental images

表1 单模态医学图像配准时浮动图像实际变换参数Tab.1 Transformation parameters of mono-modality floating images

图7 ICPFig.7 ICP

图8 RIKMCFig.8 RIKMC

表2 单模态医学图像配准浮动图像性能比较Tab.2 Performance of registering the mono-modality images

2.2 多模态医学图像配准

在该实验中参考图像选取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级。每幅参考图像变换后的浮动图像的实际变换参数如表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 多模态医学图像配准时浮动图像实际变换参数比较Tab.3 Transformation parameters of multi-modality floating images

我们首先使用ICP来完成单模态配准，然后使用RIKMC来实现单模态配准。实验结果如图13、图14、图15以及表4所示。

从表4可以看到，在多模态配准时，基于ICP的配准方法和RIKMC配准方法相比较，RIKMC配准方法整个配准时间较短；特别是图像较大时，RIKMC配准速度优势更明显，但是当图像较小，RIKMC和ICP配准速度相差不大。至于配准精度，ICP对2组MR2和CT2配准基本成功，对于1组CT1和MR1、3组 CT3和PET1、4组MR3和PET2配准配准完全失败，陷入局部最优，因此整体上讲，ICP配准失败率较高；RIKMC均能成功配准，虽然对3组CT3和PET1能成功配准，但是配准误差相对较高。因此，从整体上讲，MICP也适合多模态配准。

图13 1组配准结果Fig.13 Result figures of the first group

图14 2组配准结果Fig.14 Result figures of the second group

图15 3组配准结果Fig.15 Result figures of the third group

图16 4组配准结果Fig.16 Result figures of the fourth group

表4 多模态医学图像配准浮动图像性能比较Tab.4 Performance of registering the multi-modality images

3 结论

为了解决ICP算法配准速度慢、配准成功率低等问题，论文利用改进的K-Means聚类方法，提出了基于改进K-Means聚类医学图像配准算法。通过实验得出该算法既可用于单模态图像配准，也可以用于多模态图像配准；具有运算量少、图像配准速度较快、计算比较简单、精确度较高等特点，并且解决了图像配准容易陷入局部最优的问题。

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