数据采集与处理  2018, Vol. 33 Issue (2): 248-258   PDF    
基于边缘特征点互信息熵的医学图像配准方法
魏本征1,2 , 甘洁3 , 尹义龙4     
1. 山东中医药大学理工学院,济南,250355;
2. 山东中医药大学计算医学实验室,济南,250355;
3. 山东中医药大学第二附属医院放射科,济南,250001;
4. 山东大学计算机科学与技术学院,济南,250100
摘要: 基于互信息熵的图像配准方法已经被广泛应用于医学图像配准中,为克服互信息配准方法的不足,结合图像空间结构信息,本文提出一种基于边缘特征点互信息熵的医学图像配准方法,设计了包括互信息熵、图像空间结构和形状特征点等多信息融合的配准新测度。算法首先采用改进的形态学梯度提取医学图像边缘轮廓;然后构造了以边缘区域特征和梯度信息为基础的特征点互信息能量函数,并通过最小化能量函数来获取配准参数;最后,结合梯度下降法优化策略,实现图像配准。实验研究表明,该方法在保证了配准精度的同时,配准速度较快、鲁棒性较好、综合性能优良,具有一定的临床推广价值。
关键词: 图像配准    医学图像    互信息熵    测度函数    边缘特征    
Medical Image Registration Based on Mutual Information Entropy Combined with Edge Correlation Feature
Wei Benzheng1,2, Gan Jie3, Yin Yilong4     
1. College of Science and Technology, Shandong University of Traditional Chinese Medicine, Jinan, 250355, China;
2. Computational Medicine Lab, Shandong University of Traditional Chinese Medicine, Jinan, 250355, China;
3. Department of Radiology, Second Affiliated Hospital of Shandong University of Traditional Chinese Medicine, Jinan, 250001, China;
4. School of Computer Science and Technology, Shandong University, Jinan, 250100, China
Abstract: Image registration is a valuable technique for medical diagnosis and treatment. Due to the inferiority of image registration using maximum mutual information, a new hybrid method of multimodality medical image registration based on mutual information of spatial information is proposed. The new measure that combines mutual information, spatial information and feature characteristics, is proposed. Edge points are used as features and obtained from a morphology gradient detector. Feature characteristics like location, edge strength and orientation are taken into account to compute a joint probability distribution of corresponding edge points in two images. Mutual information based on this function is minimized to find the best alignment parameters. Finally, the translation parameters are calculated by using a gradient descent algorithm. The experimental results demonstrate the high validation precision and excellent accelerating capability of the algorithm.
Key words: image registration    medical image    mutual information entropy    measure function    edge feature    
引言

医学图像配准是医学图像处理领域中的一项重要技术,在医学基础研究和临床诊断治疗中起着越来越重要的作用。医学图像配准是使包含成像数据集中的信息最大化的一个重要步骤,其任务是对不同获取时间、不同成像设备及不同获取角度的数据进行空间上的匹配,通过匹配可以得到多源医学图像信息的空间对应关系。因此,医学图像配准算法的主要任务就是确定一个映射,使图像集在空间上关联,从而使这些图像能够在一个统一的坐标系上表示出来。

医学图像配准技术的基本框架主要包括特征空间、搜索空间、相似性测度和搜索策略4个方面。其中,作为图像配准依据的图像信息集合称之为特征空间,其主要是指从两幅待配准图像(参考图像和浮动图像)中提取具有区别力的图像信息特征,这对配准算法的运行速度和鲁棒性等性能具有重要影响。配准用图像特征通常应具有较高的普适性和稳定性,较小的特征匹配计算量,以及特征提取简单、快捷等特性。当前图像配准方法一般分为两种,一种是基于图像特征的配准方法,其特征空间主要包括点、边缘、曲线、曲面、形状或是特征描述子,如轮廓矩、Zernike矩、小波描述符及方向链码等。该类方法主要是利用如点、线段和面等图像形状共同特征构建配准目标函数,配准速度快,但特征提取困难且不稳定,误配准率较高。另一种是基于灰度的配准方法[1-2],图像的灰度值特征是其最主要的特征空间。此类配准方法主要是对图像的灰度进行操作,算法思想简单,且无需预处理,然而该配准过程计算量过大,易受到图像噪声和形状变换的影响。在基于灰度的各种配准测度中,互信息熵配准方法具有鲁棒性好、精度高等优点,但其计算量大,而且当图像受低采样、噪声、变形等因素影响时,该配准方法易出现误匹配、鲁棒性差等问题[3]

在诸多改进图像配准算法的研究中,将几何特征与像素相似性方法相结合,实现两者优势互补,从而设计出更加稳定、性能更好的相似性配准测度,是一种有效和可行的研究思路。Rangarajan等人[4]提出了一种利用互信息匹配形状特征点进行配准的策略。Pluim等人[1]通过提取图像的梯度信息,并将其视为空间信息加入互信息中,获得了较高的匹配精度。火元莲等人[5]提出一种基于轮廓特征点最大互信息的多源医学图像配准方法。卢振泰等人[6]则采用先整体后局部的刚性配准策略,提出一种基于等效子午面与互信息量的三维医学图像快速配准算法。赵海峰等人[7]提出一种基于特征点Rényi互信息的医学图像配准算法,该算法适于单模和多模医学图像配准,图像配准速度较快,鲁棒性较强。

借鉴上述研究方法,为进一步提高互信息熵类医学图像配准方法的效率和精度,基于医学图像局部特征信息,本文提出了一种基于边缘特征点互信息熵的医学图像配准方法(Feature point based mutual information method,FPMI)。

1 边缘特征点互信息熵配准方法

为达到提高医学图像配准质量和配准速度,降低图像噪声对图像配准影响的目的,基于图像轮廓、边缘特征点和互信息熵等图像局部信息,所提出的FPMI算法主要包括:计算图像熵及互信息量,提取待配准图像边缘,构建边缘特征点互信息能量函数,以及选取配准优化策略等。

1.1 计算图像熵及互信息量

熵表达的是一个系统的复杂性或者是不确定性。对于一幅图像I来说,它含有信息量的多少可以用熵表示为

$ H\left( I \right) = - \sum {{p_x}{{\log }_2}{p_x}} $ (1)

式中px表示医学图像中图像灰度强度值为I(x)的像素点在所在图像中的出现频率,即相关统计概率。

图像互信息量是指两幅图像中的共有信息量,当两幅图像实现了正确匹配的时候,它们相关的互信息量最大。所以,在选择医学图像配准用的相似性测度函数时,可以选取用两幅图像的互信息熵作为它们图像配准的依据[8]

在表示两幅医学图像XY的互信息熵MI(X, Y)时,可选用联合信息熵的形式, 即

$ MI\left( {X,Y} \right) = H\left( X \right) + H\left( Y \right) - H\left( {X,Y} \right) $ (2)

由Dobrushin公式[9]可推导得互信息熵的计算公式为

$ MI\left( {X,Y} \right) = \sum\limits_{ij} {{p_{XY}}\left( {x,y} \right){{\log }_2}\frac{{{p_{XY}}\left( {x,y} \right)}}{{{p_X}\left( x \right){p_Y}\left( y \right)}}} $ (3)

式中边缘概率分布计算公式pX(x)和pY(y)描述了两幅待配准图像XY的边缘概率分布情况。另外,待配准图像XY的联合概率分布是pXY(x, y),直接计算比较困难,可借助于两幅图像的联合直方图h(x, y)进行计算。h(x, y)表示两幅图像重叠部分图像的灰度值为(x, y)的像素对总数,其中xy分别表示两图中灰度等级,(x, y)表示两图中同一个位置的像素值,如果浮动图像溢出边界,则两幅图像中不重合的部分不参加计算,或是用图像背景灰度填充因浮动图像几何变换而出现的空白区域。设两幅图像的重叠部分像素点为(i, j) (i = 0, 1, 2, …,M-1; j= 0, 1, 2, …,N-1),其联合直方图可以表示为[10]

$ h\left( {x,y} \right) = \left[ {\begin{array}{*{20}{c}} {h\left( {0,0} \right)}&{h\left( {0,1} \right)}& \cdots &{h\left( {0,N - 1} \right)}\\ {h\left( {1,0} \right)}&{h\left( {1,1} \right)}& \cdots &{h\left( {0,N - 1} \right)}\\ \vdots&\vdots&\ddots&\vdots \\ {h\left( {M - 1,0} \right)}&{h\left( {M - 1,1} \right)}& \cdots &{h\left( {M - 1,N - 1} \right)} \end{array}} \right] $ (4)

pY(y)可以近似为两幅图像重叠部分像素对的数量(图像X中第x灰度级和图像Y中第y灰度级构成的像素对)与总的像素对数(所有的灰度的像素对数)之比。利用联合直方图,两幅图像的联合概率分布为

$ {p_{XY}}\left( {x,y} \right) = \frac{{h\left( {x,y} \right)}}{{\sum\limits_{i \in X,j \in Y} {h\left( {x,y} \right)} }} $ (5)

边缘概率分布pX(x)和pY(y)也可以通过联合概率密度求得

$ {p_X}\left( x \right) = \sum\limits_y {{p_{XY}}\left( {x,y} \right)} $ (6)
$ {p_Y}\left( y \right) = \sum\limits_x {{p_{XY}}\left( {x,y} \right)} $ (7)

通过上述公式的推导和计算,不难发现对于互信息熵的计算,可以转化成为对两幅图像的联合概率分布和边缘概率分布的计算。

假设两幅待配准图像X(i, j)和Y(i, j)的灰度级数都是256,尺寸大小为M×N,其中(i, j)是两幅图像重叠部分像素点的坐标,两幅图像中未重叠的部分与计算无关。则计算两幅图像互信息熵MI(X, Y)的算法如下。

算法1  图像互信息熵的计算算法

(1) 确定参考图像、浮动图像的坐标;初始化联合直方图hXY(x, y)数据,并令hXY(x, y)=0;

(2) 对于两幅图像的重叠区域上的每一点对(i, j),(i= 0, 1, 2, …, M-1; j= 0, 1, 2, …,N -1), 计算hXY(X(i, j), Y(i, j))=hXY(A(i, j), B(i, j))+1;

(3) 作归一化处理$p\left( {x,y} \right) = \frac{{h\left( {x,y} \right)}}{{M \times N}}$

(4) 分别计算图像X, Y的边缘概率分布${p_X}\left( x \right) = \sum\limits_{b = 0}^{255} {{p_{XY}}\left( {x, y} \right)} $${p_Y}\left( y \right) = \sum\limits_{b = 0}^{255} {{p_{XY}}\left( {x, y} \right)} $

(5) 利用式(3),可求得两幅图像的互信息熵。

1.2 图像边缘特征提取

为得到一个更加准确的医学图像轮廓,有必要消除图像噪声的干扰,从而得到尽可能完整的图像边缘轮廓。相比传统的边缘检测方法,Matherom和Serra根据几何代数及拓扑论提出的数学形态学具有更好的抑制噪声效果[11]。形态学梯度算子可表示为

$ {\rm{Grad}}\left( f \right) = \left( {f \oplus g} \right) - \left( {f \ominus g} \right) $ (8)

式中:f(x, y)为原始图像,g(x, y)为结构元素,⊕表示膨胀运算,$ \ominus $表示腐蚀运算。显而易见,形态学梯度算子能加剧输入图像的灰度级阶跃变。

为使得提取的医学图像边缘曲线更平滑,同时保证边缘特征的完整性和细节丰富性,借鉴广义形态滤波器原理,设计边缘检测算子[12]和形态学梯度滤波算子[13],如式(9, 10)所示。

$ I{\rm{Grad}}\left( f \right) = \left( {f \circ g} \right) \oplus g - \left( {f \cdot g} \right) \ominus g $ (9)
$ F{\rm{Grad}}\left( f \right) = {c_1}I{\rm{Grad}}1\left( f \right) + {c_2}I{\rm{Grad}}2\left( f \right) $ (10)

式中:c1c2的权系数大小采用最小均方自适应方法确定;IGrad1()和IGrad2()为分别与3×3的十字形和交叉形结构元素相对应的边缘检测算子。

1.3 边缘特征点互信息能量函数

由于待配准的医学图像几乎都基于共同解剖结构信息,因此在尽量不减少配准图像特征信息的前提下,为减少互信息熵的计算量,提高配准效率,基于互信息的配准策略可采用相同特征点集的互信息熵来取代整幅医学图像的互信息熵作为配准目标函数[5]

假设从两幅待配准医学图像中分别提取出的形状特征点的集合为X={Xi, i=1, 2, …, N1}和Y={Yi, i=1, 2, …, N2}, 其中XiYi是点在二维平面中的坐标位置。则基于变换参数TXY的距离测度可表示为

$ D\left( T \right) = \sum\limits_{i = 1}^N {{{\left\| {X - TY} \right\|}^q}} $ (11)

式中:N=min{N1, N2},参照K-means聚类算法[14],选取参数q =3.5。根据互信息的定义,可得点集的互信息熵

$ MI\left( P \right) = MI\left( {X,Y} \right) = \sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {{P_{ij}}\log \frac{{{P_{ij}}}}{{\sum\nolimits_{k = 1}^{{N_1}} {{P_{kj}}} \sum\nolimits_{k = 1}^{{N_2}} {{P_{il}}} }}} } $ (12)

式中:Pij为两幅图像边缘特征点XiYj的联合概率,即同时从X中选取Xi和从Y中选取Yj的联合概率,Pij=P{I=i, i∈(1, 2, …, N1), J=j, j∈(1, 2, …, N2)}。由式(12)可知,当两幅图像配准时,对应匹配点的空间位置一致,如果有NXiYj是对应匹配点,N的取值应为最大,此时图像间的相关性最大,互信息量也取得最大值。因此实际计算时,可采用如下经验计算公式[4]

$ {P_{ij}}\left( T \right) = \exp \left( { - \alpha {D_{ij}}\left( T \right) - \lambda } \right) $ (13)

式中:αλ是两个拉格朗日常数,作为辅助参数分别用于约束点匹配相似测度和概率总值。特征点互信息能量函数可表示为

$ {E_{MI}}\left( {P,T,\lambda ,\alpha ,\kappa } \right) = \alpha \left( {\sum\limits_{ij} {{P_{ij}}{D_{ij}}\left( T \right)} - d} \right) + \lambda \left( {\sum\limits_{ij} {{P_{ij}}} - 1} \right) + \\ \sum\limits_{ij} {{P_{ij}}\log {P_{ij}}} - \kappa MI\left( P \right) $ (14)

式中:噪声标准水平$d = \sum\limits_{ij} {D\left( T \right)\frac{{{\rm{exp}}\left( { - \alpha D\left( T \right)} \right)}}{{\sum\nolimits_{ij} {{\rm{exp}}\left( { - \alpha D\left( T \right)} \right)} }}} $κ为一个权重参数,且κ>0,用于调整互信息和距离测度之间的比重关系。最小化能量函数的过程如式(15)所示。

$ {\partial ^ * } = \mathop {\arg }\limits_{\left( T \right)} \min {E_{MI}}\left( {P,T,\lambda ,\alpha ,\kappa } \right) $ (15)

上述能量函数设计过程中,K-means聚类算法被用于提取医学图像特征点集,基于边缘点特征点集XY,选取聚类中心K =200,即选取200个点。之后,需找到一个最优的配准参数,使得EMI最小。

1.4 配准优化策略

在具体应用过程中,因为参考图像或浮动图像的灰度特征信息以及图像含有噪声的干扰,互信息测度函数在配准过程中存在大量局部极值点。文献[10]已经对配准过程中局部极值的产生原因和可能采取的克服办法进行了专门的研究,对此本文不再进行讨论。在多维空间,因为大量极值点的存在使得参数优化成为一个难题。从计算的角度看,在确定了图像配准目标函数和特征点集之后,图像配准就变为了一个配准参数寻优的过程。相应算法即是基于某种搜索策略和图像变换方法,在搜索空间中找最优解,求解目标函数EMI最小值的过程。最优解的求解过程就是求参考图像和浮动图像之间变换模型参数的过程[15]。一般来说,相似性测度的计算比较复杂,并且会产生大量数据,所以常采用最优化求解法来减少运算量,加快配准速度。在实际求解过程中,几乎每一种配准过程都需要最优化算法,用来搜索使目标函数取得最小值(或使用相似性测度取得最大值)的最优变换。对一个特定的配准算法,配准过程可表示为

$ {T_{{\rm{Optimal}}}} = \mathop {\arg }\limits_{\left( T \right)} \min f\left( {T\left( {{X_{\rm{F}}}} \right),{X_{\rm{R}}}} \right) $ (16)

式中:T为配准变换,XR为参考图像,XF为浮动图像,f为需要找到最小值的目标函数。

图像配准参数的优化主要有两个要求:全局寻优和快的优化速度。目前常用的优化方法主要有:梯度下降法、Powell算法和下降单纯形法,以及模拟退火算法、遗传算法和粒子群优化算法等。其中前3种算法的优点是收敛速度快、计算量小,缺点是易收敛到局部最优;后3种算法的优点是全局寻优能力强,缺点是收敛速度慢、计算量大[16]

本文选用梯度下降法作为配准测度函数的优化策略,其优化算法求解过程和执行步骤如下。

基于梯度的最优化方法通常用来确定搜索方向,在该方向目标函数的值局部递减。设xn维向量x=[x1, x2, …, xn]T,目标函数f(x)的梯度向量能够表示为

$ \nabla f\left( \mathit{\boldsymbol{x}} \right) \equiv g\left( \mathit{\boldsymbol{x}} \right){\left[ {\frac{{\partial f}}{{\partial {x_1}}},\frac{{\partial f}}{{\partial {x_2}}}, \cdots ,\frac{{\partial f}}{{\partial {x_n}}}} \right]^{\rm{T}}} $ (17)

其二阶偏导数能够用一个Hessian矩阵表示

$ {\nabla ^2}f\left( \mathit{\boldsymbol{x}} \right) \equiv H\left( \mathit{\boldsymbol{x}} \right){\left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}f}}{{{\partial ^2}{x_1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_n}}}}\\ \vdots&\ddots&\vdots \\ {\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_n}}}}& \cdots &{\frac{{{\partial ^2}f}}{{{\partial ^2}{x_n}}}} \end{array}} \right]^{\rm{T}}} $ (18)

函数f(x)以xk的泰勒展开近似表示为

$ f\left( {{\mathit{\boldsymbol{x}}_k} + \mathit{\boldsymbol{x}}} \right) \approx f\left( {{\mathit{\boldsymbol{x}}_k}} \right) + {\left( {\nabla f\left( {{\mathit{\boldsymbol{x}}_k}} \right)} \right)^{\rm{T}}}\mathit{\boldsymbol{x}} + \frac{1}{2}{\mathit{\boldsymbol{x}}^{\rm{T}}}{\nabla ^2}f\left( {{\mathit{\boldsymbol{x}}_k}} \right)\mathit{\boldsymbol{x}} = f\left( {{\mathit{\boldsymbol{x}}_k}} \right) + g{\left( {{\mathit{\boldsymbol{x}}_k}} \right)^{\rm{T}}}\mathit{\boldsymbol{x}} \\ + \frac{1}{2}{\mathit{\boldsymbol{x}}^{\rm{T}}}H\left( {{\mathit{\boldsymbol{x}}_k}} \right)\mathit{\boldsymbol{x}} $ (19)

算法2  基于梯度的最优化算法

(1) 初始化。设置迭代器k=0,并初始化向量x,计算f(xk)。

(2) 验证收敛标准。如果满足了收敛标准,则停止最优化过程,得到解xk

(3)计算搜索方向。计算搜索方向向量Pk

(4) 计算步长lk。确定满足f(xk+lkPk)<f(xk)的正标量lk的值。

(5) 更新变量。令xk+1= xk+lkPk,计算f(xk+1),令k=k+1,返回步骤(2)。

搜索方向Pk的计算和步长lk的寻找是基于梯度优化算法的两个主要问题,用不同的方法计算搜索方向就产生了各种不同的基于梯度的方法。

2 图像配准过程

Rangarajan算法[4]在定义了基于最大互信息熵后,又以其为基础,通过引入辅助变量作代数变换定义了一个能量函数EMI。两幅图像配准时分别作为参考图像和浮动图像。通过对配准图像的空间变换,使得两幅图像的像素点尽可能对应匹配,此时两幅图像的互信息达到最大值,且EMI取得最小值。在计算过程中,一般选取边缘轮廓清晰、图像特征提取稳定的图像作为参考图像,另一幅作为浮动图像。

结合两幅图像边缘特征点互信息熵的计算方法,采用梯度下降法作为图像配准优化算法。算法流程如下。

算法3  FPMI图像配准算法

(1) 用形态学梯度滤波算法获取两幅图像轮廓,并采用K-means聚类算法提取图像边缘特征点, 得到形状特征点集XYX为参考图像点集, Y为浮动图像点集;

(2) 初始化,Ρij=(N1N2)-1, α=αminλ=λ0κ=κminT=0;

(3) 确定算法的搜索空间,ΔT=(Δx, Δy, Δθ),Δα,Δκ

(4) 初始化由ΔT,Δα和Δκ构成梯度优化算法的5维向量x;

(5) 计算待配准图像的梯度损失函数值(两幅图像的EMI);

(6) 根据式(15, 19)更新x

(7) 迭代结束条件通常为全局最优解小于最小允许误差, 或达到一个预设最大代数Gmax,如未达到该结束条件则返回步骤(5);

(8) 达到迭代结束条件,配准完成。

3 实验结果与分析

本文将FPMI算法与基于最大互信息熵等算法在配准精度、运行时间等算法性能方面进行比较。用Matlab 7.0在PC机上(P4 2.8 GHz CPU, 1 GB内存)分别对多源医学图像进行配准实验。利用直方图估计概率密度分布时采用64位灰度级。配准成功率设定为配准成功的次数与实验总数的比值,且每套数据独立进行100次实验。利用形态学梯度滤波算法提取边缘时采用的结构元素为{0 1 0; 1 1 1;0 1 0}T和{1 0 1;0 -1 0; 1 0 1}T。优化算法均采用梯度下降法,最大迭代步数200次,算法初始化的搜索空间为:x方向的搜索范围x∈[-20, 20],y方向的搜索空间y∈[-20, 20],旋转的搜索空间角度θ∈[-15°, 15°]。

3.1 模拟配准测试实验 3.1.1 人工合成图像配准实验

为检验配准算法的正确性,本文采用人工合成图像及其变换图像作为配准实验的参考和浮动图像,对所提FPMI算法做正确性验证配准实验。待配准图像如图 1所示, 图 1(a)为参考图像,图 1(b)为浮动图像,两幅图像中白色矩形区域有一定的旋转角度和平移量偏差。其变换参数(Δx, Δy, Δθ)为(40, 20, 20°)。

图 1 人造合成配准测试实验图像 Fig. 1 Synthetic images for registration test

实际计算过程中,配准成功的标准设定为:如果计算出的变换参数(Δx, Δy, Δθ)和平移量与真实值相差1个像素,且旋转角度偏差值在1°范围内,即认为配准成功。100次算法验证实验,成功率为100%。配准结果如表 1所示,FPMI算法平移配准精度在0.2像素以内,旋转角度精度在0.1°以内。实验结果表明了FPMI算法的正确性。

表 1 人工合成图像FPMI配准实验结果统计表 Tab. 1 Registration results of FPMI on synthetic image

3.1.2 CT-MRI配准实验

为检验FPMI算法的有效性,本文对CT-MRI图像进行模拟配准对比实验,结果如图 2所示,图 2(a)为CT图像,作为参考图像,图 2(b)为MRI图像,作为浮动图像。每幅图像的大小为512像素×512像素,图像取自BIS医学图像库(http://nova.nlm.nih.gov/Mayo)。在实验中,图像配准分别采用最大互信息熵算法(MI算法)、具有代表性的特征互信息量配准算法(Rangarajan算法[4])以及本文的FPMI算法。各算法均独立运行100次,实验结果如图 2(c, d)表 23所示。其中,图 2(c)为FPMI算法求解的MRI变换后图像,图 2(d)为配准后的融合图像;表 2是算法模拟配准变换参数(Δx, Δy, Δθ)的均值统计表;表 3是3种算法的配准指标统计表,包括4个指标:平均互信息量、互信息量标准差δ、单次平均计算时间和配准成功率。

图 2 脑部MRI-CT图像模拟配准示意图 Fig. 2 Simulation registration of brain MRI and CT image

表 2 算法模拟配准变换参数均值统计表 Tab. 2 Simulation registration transformation parameter

表 3 CT-MRI模拟配准各算法指标统计表 Tab. 3 Statistics of simulation registration algorithm on CT-MRI

从实验结果可以看出,对模拟CT、MRI图像进行配准时,FPMI的平均配准偏差较小,配准参数(Δx, Δy, Δθ)的精度高于另外两种方法;FPMI算法的平均配准互信息量略高于MI算法的平均互信息量,这在一定程度上保证了算法配准的准确度。MI算法成功率略低,为98%;Rangarajan算法和FPMI两种算法成功率均可以达到100%,但FPMI算法运算速度快,其运算速度约为Rangarajan算法的1.5倍,MI算法的2.4倍。因此,FPMI算法的综合性能优于另外两种算法。

3.2 测度函数曲线锐度测试实验

为对FPMI配准测度函数进行性能测试,本文选取MRI脑图像进行实验。实验所用的两幅MRI脑图像取自BrainWeb(http://www.bic.mni.mcgill.ca/brainweb/.),分别是大小为181×217×181的T1加权和PD加权的脑部MRI图像,其体素大小为1.0mm3图 3为测试实验图像及边缘特征点集示意图。如图 3(a)所示,将MRI-T1图像沿水平方向上平移一定像素后与MRI-PD图像(图 3(c))之间计算测度值,观察峰值的变化,像素平移的范围是[-10, 10],负号表示向左平移,正号表示向右平移;同样将MRI-T1图像绕中心旋转一定的角度后与MRI-PD图像计算测度值,再观察峰值的变化,角度旋转的范围是[-100°, 100°],负号表示顺时针旋转,正号表示逆时针旋转。图 3(b, d)分别为图 3(a, c)的图像边缘特征点集示意图。

图 3 测试实验图像及边缘特征点集示意图 Fig. 3 Test images and the edge feature point set

本文选用的配准测度函数是互信息量和FPMI的边缘特征点互信息量,结果如图 4所示。图 4(a, b)分别表示两测度函数值随图像在水平方向平移像素和角度旋转时的变化曲线,其中,横轴分别表示水平方向像素的平移量和旋转角度,纵轴均表示测度函数值。从图 4可得,两图像在水平平移为0像素及旋转角度为0°时测量函数均有最大值,FPMI的配准方法中测度函数的曲线峰值比MI更加尖锐,更有利于最优变换参数的选取。

图 4 配准测度函数变化示意图 Fig. 4 The comparison chart of registration measure function

3.3 多模医学图像配准实验

为了检验FPMI新算法在多模图像配准上的配准精度,本文选取两组同一病人的脑部图像进行配准实验,如图 5所示。每幅图像的大小均为256×256,并且选取MRI图像(图 5(ac))作为参考图像,CT(图 5(b))和SPECT(图 5(d))分别作为浮动图像,将本文算法与MI算法和Rangarajan算法进行了对比,3种算法独立测试50次。表 4给出了MRI-CT和MRI-SPECT的两组图像的配准结果。

图 5 多模医学配准图像 Fig. 5 Multi-modal medical image registration

表 4 两组多模医学图像配准结果统计表 Tab. 4 Multi-modal medical image registration result

在MRI-CT的配准过程中,FPMI算法和Rangarajan算法均出现一次误配,成功率达到了98%,而MI的成功率是94%;在运算效率方面,FPMI的运算效率约是MI算法的3.5倍,是Rangarajan算法的1.8倍。

在MRI-SPECT的配准实验中,3种算法均出现了不同程度的误配准现象,FPMI算法成功率达到了90%,而MI算法只取得了72%的成功率。这是由于SPECT图像灰度分布较为模糊,且其形状轮廓与MRI图像差别较大,因此导致了分别基于灰度、轮廓及边缘特征点的3种配准算法的准确率出现下降。在运算效率方面,FPMI的运算时间达到了1 027 s,其效率约为MI的2.3倍,Rangarajan算法的1.9倍。与CT-MRI配准实验相比,在组织结构更为复杂的医学图像配准中MI和Rangarajan算法的运算效率降低幅度更大,误配率更高。

因此,对比实验结果可知,相对于另外两种算法,FPMI算法在配准精度和运算效率上具有较为明显的优势。且将其应用于边缘轮廓更为清晰的医学图像配准时,其算法的性能优势会发挥得更好。

从上述实验可以看出,与MI算法及Rangarajan算法相比较,FPMI算法在配准精度和运算效率上都有一定优势。实验结果还表明,FPMI算法中采用MGZ算法提取图像轮廓特征较为稳定和准确,所以此方法对于图像轮廓相对清晰的医学图像是一种较理想的配准方法,但该方法对数据的缺失及图像轮廓的形状变异较为敏感。此外,相对于同类算法,在配准过程中本算法无需执行粗配准和精配准两个过程,因此算法的配准过程简单、实用性更强。

4 结束语

基于互信息熵配准测度函数,本文设计了一种基于多特征信息融合的边缘特征点互信息配准测度,并利用梯度下降法优化配准过程。所提出的FPMI算法能反映医学图像的空间形状特征,可大幅减少互信息的计算量,大大缩短配准时间,并可有效克服局部极值问题。此外,本文还研究用梯度下降法进行医学图像配准的参数寻优,以进一步提高配准的效率和准确度。研究结果表明,该方法能够准确、快速地处理刚性配准问题,在多源医学图像的配准中具有一定的先进性,是一种稳健的自动配准方法,可用于临床研究,具有较高的实用价值。但是,本文算法并不能适用于所有各类医学图像,因此下一步工作将针对不同医学图像的配准需求研究更具有针对性的配准测度算法。

参考文献
[1]
Pluim J P W, Antoine M J B, Viergever M A. Image registration by maximization of combined mutual information and gradient information[J]. IEEE Transactions on Medical Imaging, 2002, 19(8): 809-814.
[2]
陈庆芳, 吴小俊. 基于分块互信息和量子粒子群算法的图像配准[J]. 数据采集与处理, 2011, 26(4): 473-477.
Chen Qingfang, Wu Xiaojun. Image registration based on block mutual information and quantum-behaved particle swarm optimization[J]. Journal of Data Acquisition & Processing, 2011, 26(4): 474-477.
[3]
何彦杰, 周焰, 刘超, 等. 结合PCM和MIE的多模态图像配准方法[J]. 空军预警学院学报, 2010, 24(5): 372-375.
He Yanjie, Zhou Yan, Liu Chao, et al. Method of multimode image registration based on PCM and MIE[J]. Journal of Air Force Radar Academy, 2010, 24(5): 372-375.
[4]
Rangarajan A, Chui H, Duncan J S. Rigid point feature registation using mutual information[J]. Med Image Anal, 1999, 3(4): 425-440. DOI:10.1016/S1361-8415(99)80034-6
[5]
火元莲, 齐永锋, 宋海声. 基于轮廓特征点最大互信息的多源医学图像配准[J]. 激光与红外, 2008, 38(1): 96-98.
Huo Yuanlian, Qi Yongfeng, Song Haisheng. Multi-modality medical image registration based on mutual information of feature points[J]. Laser & Infrared, 2008, 38(1): 96-98.
[6]
卢振泰, 冯衍秋, 冯前进, 等. 基于等效子午面与互信息量的医学图像配准[J]. 计算机学报, 2009, 32(8): 1611-1617.
Lu Zhentai, Feng Yanqiu, Feng Qianjin, et al. Medical image registration using equivalent meridian plane and mutual information[J]. Journal of Computer, 2009, 32(8): 1611-1617.
[7]
赵海峰, 陆明, 卜令斌, 等. 基于特征点Rényi互信息的医学图像配准[J]. 计算机学报, 2015, 38(6): 1212-1220.
Zhao Haifeng, Lu Ming, Bu Lingbin, et al. Medical image registration based on feature points and Rényi mutual information[J]. Chinese Journal of Computers, 2015, 38(6): 1212-1221. DOI:10.11897/SP.J.1016.2015.01212
[8]
Penney G P, Weese J, Little J A, et al. A comparison of similarity measures for use in 2D-3D medical image registration[J]. IEEE Trans Med Imag, 1999, 17(4): 586-595.
[9]
Wells W M, Viola P, Atsumi H, et al. Multi-modal volume registration by maximization of mutual information[J]. Med Image Anal, 1996, 1(1): 35-51. DOI:10.1016/S1361-8415(01)80004-9
[10]
Maes F, Collignon A, Vandermeulen D, et al. Multimodality image registration by maximization of mutual information[J]. IEEE Trans Med Imag, 1997, 16(2): 187-198. DOI:10.1109/42.563664
[11]
Serra J. Image analysis and mathematical morphology[M]. New York, NY, USA: Academic Press, 1982: 1-610.
[12]
魏本征, 赵志敏, 宋一中. 基于IPSO和综合信息的医学图像配准新方法[J]. 光电子·激光, 2009, 20(9): 1271-1274.
Wei Benzheng, Zhao Zhimin, Song Yizhong. A novel algorithm for multimodality medical image registration based on IPSO algorithm and hybrid information[J]. Journal of Optoelectronics Laser, 2009, 20(9): 1271-1274.
[13]
魏本征, 赵志敏, 华晋. 基于形态学梯度和Zernike矩的亚像素边缘检测方法[J]. 仪器仪表学报, 2010, 31(4): 838-844.
Wei Benzheng, Zhao Zhimin, Hua Jin. Sub-pixel edge detection method based on improved morphological gradient and Zernike moment[J]. Chinese Journal of Scientific Instrument, 2010, 31(4): 838-844.
[14]
Hamerly G, Elkan C. Alternatives to the k-means algorithm that find better clusterings[C]//Proc of the ACM Conference on Information and Knowledge Management, CIKM-2002. [S. l. ]: CIKM, 2002: 600-607.
[15]
王伟. 医学图像非刚性配准算法研究[D]. 大连: 大连理工大学, 2012.
Wang Wei. Research on non-rigid registration algorithm of medical images[D]. Dalian: Dalian University of Technology, 2012. http://cdmd.cnki.com.cn/Article/CDMD-10141-1012393375.htm
[16]
龚晓彦. 基于互信息的医学图像配准算法研究[D]. 秦皇岛: 燕山大学, 2010.
Gong Xiaoyan. Research on medical image registration based on mutual information[D]. Qihuangdao: Yanshan University, 2010. http://cdmd.cnki.com.cn/Article/CDMD-10422-2008191406.htm
基于边缘特征点互信息熵的医学图像配准方法
魏本征 , 甘洁 , 尹义龙