基于MRF模型的多模态图像配准技术研究
袁桂霞
摘 要： 在医学图像处理领域中，医学图像配准技术极其重要，其价值体现在临床医学中对图像处理技术的应用。在解决多模态图像配准的相关问题时，基于互信息方法的应用最广泛，但在某些特定的应用中该方法受到的约束仍然较多。针对这一情况，提出一种新的医学图像配准算法，模态变换的引入作为此研究算法的基础，之后新的马尔可夫能量函数则根据两幅通过模态变换后的图像矩阵以及原配准图像得以构建。同时，为了优化能量函数引入了一种改进的梯度下降算法，从而得到配准结果。最后，运用不同的医学图像进行配准实验来验证该算法，通过实验证明该配准算法具有良好的有效性及抗噪性能。
关键词： 图像配准； 模态变换； 马尔可夫随机场； 梯度下降算法； 非刚体配准； 图像矩阵
中图分类号： TN911.73?34； TP311.17 文献标识码： A 文章编号： 1004?373X（2018）01?0057?05
Abstract： In the field of medical image processing， the medical image registration technology is extremely important. Its value is embodied in the application of image processing technology in clinical medicine. For solving the relevant issues of multimodal image registration， the method based on mutual information is most?widely used， but constrained in some specific applications. In view of this situation， a new medical image registration algorithm is proposed. The modal transformation is introduced as the basis of the studied algorithm. After that， the new Markov energy function is constructed according to the two image matrixes after modal transformation and original registration image. In order to optimize the energy function， an improved gradient descent algorithm is introduced to get the registration results. The algorithm was verified with registration experiments of different medical images. The experimental results show that the registration algorithm has perfect effectiveness and high anti?noise performance.
Keywords： image registration； modal transformation； Markov random field； gradient descent algorithm； non?rigid registration； image matrix
0 引 言
图像配准是将图像信息进行汇总并对齐图像的过程，多模态图像配准即为一种将多模态图像进行空间几何关系对齐的技术。近年来，将临床医学和图像处理技术相结合，更加促进了医学图像配准技术的发展。临床上的大多数需求可以被目前已有的算法满足，但是考虑到图像个体之间内容的差异性比较大，而且包含的信息也比较多，所以现今的各种配准算法中都有一定的限制因素存在，一般只能在特定情况下满足配准。考虑到在临床医学方面图像配准具有广泛的应用，而每一种不同的应用均需要不同的形变场，因此关于医学图像配准的问题并没有一种具有普适性的算法解决。针对这一情况，本文提出了一种新的医学图像配准算法。
1 MRF随机场的组成
在图像分析中，马尔可夫随机场的问题可以视为一个建模问题，具体到图像配准中即构造一个有待优化的能量函数问题。基于马尔可夫随机场构造的能量函数（或称为目标函数）不仅和其相邻的像素点有关，而且还需要像素点之间的其他信息，这些信息可以视为一个约束项。一个能量函数可以从以下步骤推导得出。
定义[I（x）]和[J（x）]是二维空间中两幅大小相同的图像，[T（·）]是每一个像素点[x]从固定图像[I（x）]到[J（x）]的一个变换模型，对于两幅给定的图像来说，定义一个观测模型如下所示：
[I（x）=JT（x）+n（x）] （1）
式中：[n（x）]是一个高斯白噪声模型，其对应的噪声模型数学期望为[x～（0，σ2）]，[σ]是噪声模型的标准差。
将两幅图像和变换模型视为随机变量，[s]和[t]表示像素位置，则可以通过图[?=（v，ε）]表示马尔可夫随机场模型中离散标记的问题。定义[?=（1，2，…，L）]为图像的一个标号集合，则所有的标号集合可以定义为：
[X=xss∈V] （2）
式中：[xs]是[?]中的一个元素。给定一个包含平移、翻转、旋转和在[x，y]方向上的尺度变换的向量[θ，]可以得到标号集合的马尔可夫随机场模型的能量函数：
[Emrf（Xθ）=Edata+Esmooth] （3）
1.1 数据项
式（3）反映了两个相对应的像素点及其邻域像素点之间的关系。通常来说，在整个定义域[Ω]中，构造一个解决配准问题的模型是定义一个数据项。这个数据项能够表示在固定图像和浮动图像之间的距离函数。根据式（1）和式（2），图像[I（x）]的噪声模型的期望是[JT（x）]，方差是[σ2，]因此马尔可夫随机场模型的数据项可以表示为：
[Edata=s∈Vθs（xs）=12σ2s∈VI（xs）-JT（x）2] （4）
参数[σ]是作为测度项的调节参数引入的，它的具体值取决于人为选定的噪声模型以及实际的应用。尽管只采用数据项进行优化的方法在某些特定应用中是可行的，但是这并非一种十分通用的方法。这一方法的局限性体现在：某些应用中方程的解不存在；某些情况下对该问题的解不惟一；某些问题中，方程的解对原始数据的依赖不具有连续性。
1.2 平滑项
上述提到的所有限制都可以归结成一个不适定问题（ill?posed problem），这一问题最直接的解决方法是采用一个平滑项做约束。对于低层级的马尔可夫随机场问题，这一项可以视为一个归一化项（regularization）。它假定了在相同的邻域系统中相邻像素点的物理特性是不会突然发生变化的。因此，对于一个邻域系统[Ni]而言，其平滑项为：
[Esmooth=λs，t∈Niθstxs，xt=λs∈Vt∈Niθstxs，xt] （5）
式中：[λ]表示一个权重值，其值通常根据不同的配准形式以及配准所采用的不同函数[θst]来选取不同的常数值。函数[θst]的选择有很多种，比如分段函数、多项式函数、对数函数等。鉴于这些函数在某些情况下会加大计算的复杂度，本文采用一种线性函数来简化优化过程。采用一种间断自适应（Discontinuity Adaptive， DA）模型来构造函数的平滑项，而且函数[θ（x）]需要满足以下方程式：
[limx→∞θ（x）=limx→∞2xf（x）=C] （6）
式中：[C]是一个正整数；函数[f（x）]是一个在定义域中的偶函数。一般来说，马尔可夫随机场模型通常是指二阶的随机场模型。对于一个给定的概率模型，马尔可夫随机场问题可以看成是一个解决最大后验概率的问题。在上述讨论条件不变的情况下将固定图像、浮动图像和变换模型视为随机变量，则根据三者之间的统计相关性可得：
[pI（xs）T，J（xt）=pIJT（xt）] （7）
但是对于一个随机变量[X，]上述的分布往往难以求得，此时，可以考虑运用Hammersley?Clifford 定理将马尔可夫随机场模型转换为吉布斯分布：
[p（Xθ）∝e-E（Xθ）] （8）
根据式（8）可得，将马尔可夫的最大后验概率的相关约束问题换一种思路进行考虑，即用最小化能量函数[minXE（Xθ）]的相关约束问题来代替则更容易解决。
2 需求分析和系统设计
为了更大限度地提高医学图像配准的精度，本文在对目标函数进行优化之前引入了一个基于直方图计算的模态变换方法。
假设有两幅不同模态的图像[I（x）]和[J（x）]已经归一化，同时[I（x）]，[J（x）]∈[0，1]，[x]为像素点位置，通过对每一个像素点进行一个遍历，每遍历一个点就更新一次联合直方图，直到最后一个点遍历完，得到图像[I]和[J]的联合直方图[H（I，J）]。
[HI（x）N，J（x）N=HI（x）N，J（x）N+1] （9）
式中：[N]是联合直方图的样条数；[x]是当前像素点的位置；[c]表示小于[c]的最大正整数值；联合直方图[H（i，j）]是一个二维矩阵；灰度值[i]在图[I]以及灰度值[j]在图像[J]中相关点具体的数目用[H（i，j）]来表示。
运用式（10）和式（11）能够将图像[I]和图像[J]的灰度值表示形式的图像矩阵转换得大致相似。通过迭代过程，搜索图像[I]在图像[J]中重叠次数最多的一个像素点。在遍历了所有像素点之后，得到一个新的图像矩阵，也即经过模态变换后的图像矩阵。同理，可以由图像[J]得到[JT。]
[IT（x）=argmaxjHI（x）N，jN] （10）
[JT（x）=argmaxiHiNN，J（x）N] （11）
图1给出了利用MRI?T1序列脑部图像和CT图像进行模态变换的结果图。图1c）和图1d）是进行模态变换前后的联合直方图，从图中可以明显看到在进行模态变换前，联合直方图上的点分布比较松散，变换之后的点分布紧凑，且向图中的斜对角线靠拢，这表明模态变换的方法减少了图像簇的分散性，在某种程度上可以提高图像的配准精度。
3 基于MRF的模态变换多模态图像配准
3.1 MRF能量函数
考虑到图像[I]以及图像[J，]在图像的定义域[Ω=1，m×1，n]已经给定的情况下，将得到的4个图像矩阵结合可以得到如下数据项：
[Edata=12σ21s∈VIT（xs）-JT（x）2+12σ22s∈VI（xs）-JTT（x）2] （12）
式中：[σ1]和[σ2]是高斯白噪声模型的标准差。在实验中，为了简化计算过程，取[σ1=σ2=σ]。通常来说，数据项来自于对距离函数的定义，它通常处理的是大多数的线性问题。为了提高配准精度以及配准过程的稳定性，本文引入了另外一项表达式表示数据项，更多的细节信息的获得來自于之前转换后的图像矩阵[IT]和[JT。]
能量函数只和局部形变有关，但是在优化过程中仍然有可能会产生局部极值的情况，因此本文选择了一个基本函数作为平滑项的势能函数：
[θ（x）=ln（1+x2）] （13）
[f（x）=11+x2] （14）
式中：函数[θ（x）]和[f（x）]均满足表达式（6）的条件。由于在形变过程中，平滑项更多的是和图像像素的邻域系统有关，同时结合表达式（5）、（12）和（13），得到平滑项的相关计算表达式如下：
[Esmooth=λs∈Vt∈Niln1+xs-xt2] （15）
平滑项的引入能够处理非线性问题，同时平滑形变场，从而可以在优化过程中减少极值对配准结果的影响。
在得到数据项和平滑项之后，根据马尔可夫随机场模型理论可以得到一个新的能量函数模型：
[Emrf=Edata+Esmooth=12σ21s∈VIT（xs）-JT（x）2+ 12σ22s∈VI（xs）-JTT（x）2+λs∈Vt∈Niln1+xs-xt2] （16）
值得注意的一点是，当图像完全配准之后，马尔可夫随机场的能量函数值[Emrf]将会最小。通常来说，一个马尔可夫随机场模型指的是一个二阶模型，在计算过程中需要求像素点相邻的八个像素点（一阶系统中，斜对角线上的像素点不考虑在内）。由于二阶邻域系统计算了更多的图像信息，从而可以得到更精确的配准结果，本文中的实验均采用二阶邻域系统。
3.2 优化算法
对于医学图像配准来说，优化方法主要分为两大类：参数化和非参数化配准。因为非参数化配准广泛应用于非刚体医学图像配准中，本文主要讨论这一种优化方法。优化方法的通用模型如下：
[xi+1=xi+ai?di，i=0，1，2，…] （17）
式中：[xi]和[xi+1]分别为第[i]次和第[i+1]次迭代过程中的优化值；[di]表示当前优化过程的搜索方向，搜索方向既可以为正值也可以为负值。沿着搜索方向，[ai]是一个在每次迭代过程中用来控制步长的标量尺度系数。
一个最基本的优化方法是梯度下降方法，它是众多无约束最优化方法中的一种。这种方法迭代简单，而且相对应的配准结果也较为准确。但是这种算法在优化时并不是逐次向最优项靠近，而且由于迭代步长较小，所以迭代时间相对较长。因此，通过用图形的子集计算互信息的差分方法来提高多模态图像配准的计算时间。这一方法用预定义的一个衰减函数取代了常数值作为每次迭代的增益。其中衰减函数的表达式如下：
[ai=ai+Aα] （18）
式中：[a>0，A≥1，0<α≤1。]这一函数能够更好地选择步长值，同时也减少了优化过程中的迭代次数。在此基础上，用自适应随机梯度算法（Adaptive Stochastic Gradient Descent，ASGD）对该算法进行改进。这一方法通过计算[di+1]和[di]内积值来调整步长[ai]的大小：
[xi+1=xi-ai?di，i=0，1，2，…] （19）
[ai+1=ai?f-dTi?di-1+] （20）
式中：[[x]+]表示[x]和0进行比较后取较大的数，即[max（x，0）；]非线性sigmoid函数用函数[f]进行表示，该函数的一般表达式为[f（x）=11+e-xw]。在迭代过程中，若目标函数的当前值小于最后一次迭代值，则保持当前步长不变，否则当前步长将被式（20）的步长值取代。
通过将上述两种优化方法进行结合的方法改进优化算法。在同一个配准实验中，这一方法提供了一个更加灵活的选择步长的方式，同时使得在优化过程中减少了迭代次数，且配准输出的结果和文献中的算法具有可比性。算法中的步长值[ai+1]定义为：
[ai+1=ai?f-dTi?di-1+，Δ≥εα?ai，Δ<ε] （21）
式中：[α]是一个在定义域为[0.9，1]的常数；[Δ]是[ai]和[ai+1]差值的绝对值[ai-ai+1]；[ε]是一个根据经验而选取的值，通常它是一个小于1的常数。将梯度[di]由[ΔEmrf]代替，即得到最终的优化算法：
[xi+1=xi-ai?ΔEmrf（xi）] （22）
由式（22）可以看到，对于能量函数以及能量函数的导数，通常需要逐点进行求解，这样对于较大尺寸的图像进行配准时，会相当耗时。因此，文中所做的实验均采用一种多分辨率（multiresolution）的配准方法。这一方法的引入可以避免配准过程中遇到的局部极值，从而提高配准精度。将这一方法应用到研究中，首先对一个16×16的子网格像素点进行配准，然后将得到的形变场和原始图像调整到32×32的子网格再进行配准，以此类推，直到配准过程中这一网格达到待配准图像的原始大小，配准结束。图2中，三条不同的曲线分别表示在每一个层级优化过程的迭代情况。
为了验证本文优化方法的有效性，本文将上述方法和前人所用的研究方法进行了对比。实验中除了优化方法以外，所有的參数都设定为相同值。在式（18）和式（21）中，参数[A=10，a=1，α=0.6，ε=0.1。]同时，多分辨率方法设定层次为三层。实验用脑部图像参见图1。迭代过程中两种优化方法在迭代过程中的曲线如图2所示。从图2a）中可以看出，Myronenko方法在前20～40次迭代时陷入了局部极值，且这一情况在优化的最后一个层级仍然存在。而本文提出的优化算法不仅在迭代次数上少于Myronenko方法（本文算法的迭代次数少于100次），而且在整个迭代过程中曲线较为平滑，较大程度上解决了图2a）出现的局部极值问题。
4 算法验证
为了验证本文所提算法的可行性以及配准效率，本研究做了一系列图像配准实验。当配准过程开始时，需要应用一个全局性的仿射变换以增加配准的准确性。初始的仿射变换参数[x=][0，0，0，100，100，0，0]分别代表在[x，y]方向上的平移、旋转、尺度变换以及裁剪（shearing）。之后，一个三层的多分辨率方法同时应用于固定图像和浮动图像，所有图像的灰度值均归一化为[0，1]的区间。在每一个层级的迭代过程中，都应用模态变化的方法。之后，优化方法对已经构建好的能量函数进行优化，能量函数的参数设定为[σ1=σ2=σ=1]。所有最大迭代次数都设定为200次。
在第一个刚体实验中，选取了MRI?T1序列脑部图像和MRI?PD序列脑部图像分别作为固定图像和浮动图像，可参见图1a）和图1b）。浮动图像是在固定图像的基础上人为旋转了10°，同时在[x]和[y]方向上分别有13 mm和17 mm的位移。图像大小均为221×257。图3给出了两幅实验图像的配准后图像以及配准前后的棋盘式叠加图。通过棋盘图可以看到图像边缘有了很好的配准效果。
实验2中，验证了3种不同的非刚体形变情况下，本文算法的配准效果。如图4所示为本实验的配准用图及配准结果。图中的每一列分别为固定图像、浮动图像、配准后图像以及配准前后的棋盘格叠加图。从配准前后叠加图的效果来看，本文提出的算法输出了准确的配准结果，图中从上到下的图像模态分别为MRI?T2/CT、MRIT2/MRIT1和CT/MRI?T2。
5 结 语
本文在结合MRF随机场和模态变换的基础上，构造了一个新的待优化能量函数。此能量函数由数据项和平滑项构成。数据项包含了4个图像矩阵，分别为固定图像和浮动图像，以及由这两幅图像进行模态变换后得到的2个图像矩阵。增加两幅图像矩阵可以增加图像的细节，通过获得更多的图像信息从而提高图像配准的精度，平滑项选取了自然对数函数。
在优化方法的选取上，本文改进了梯度下降算法，将两种优化方法进行结合，使优化过程中的迭代步长选择更加灵活。通过实验验证，此方法能够有效地减少图像迭代次数，大大降低了优化过程陷入局部极值的情况。为了验证本文算法的配准精度以及抗噪性能，本文用MRI不同序列脑部图像和CT图像做了一系列的配准实验。实验结果表明，本文提出的配准方法配准效果更为准确，同时抗噪性能更好。
参考文献
[1] 张冉，王雷，夏威，等.2D/3D图像配准中的相似性测度和优化算法[J].激光与红外，2014（1）：98?102.
ZHANG Ran， WANG Lei， XIA Wei， et al. Comparison of similarity measurement and optimization [J]. Laser & infrared， 2014（1）： 98?102.
[2] 王雷，黄晨雪，王高波.改进的MRF彩色图像分割算法[J].计算机工程与应用，2016，52（13）：191?194.
WANG Lei， HUANG Chenxue， WANG Gaobo. Improved MRF color image segmentation algorithm [J]. Computer engineering and applications， 2016， 52（13）： 191?194.
[3] 张静亚，王加俊.一种改进的非刚性医学图像配准算法[J].计算机应用研究，2015，32（4）：1261?1264.
ZHANG Jingya， WANG Jiajun. Improved non?rigid medical image registration algorithm [J]. Application research of computers， 2015， 32（4）： 1261?1264.
[4] 金斌，周伟，丛瑜，等.一种基于局部不变特征的SAR图像配准新算法[J].哈尔滨工业大学学报，2014，46（11）：112?118.
JIN Bin， ZHOU Wei， CONG Yu， et al. A novel and efficient algorithm using local invariant feature for SAR image registration [J]. Journal of Harbin Institute of Technology， 2014， 46（11）： 112?118.
[5] 赵雪梅，李玉，赵泉华.结合高斯回归模型和隐马尔可夫随机场的模糊聚类图像分割[J].电子与信息学报，2014，36（11）：2730?2736.
ZHAO Xuemei， LI Yu， ZHAO Quanhua. Image segmentation by fuzzy clustering algorithm combining hidden Markov random field and Gaussian regression model [J]. Journal of electronics & information technology， 2014， 36（11）： 2730?2736.
[6] 徐胜军，韩九强，何波，等.融合边缘特征的马尔可夫随机场模型及分割算法[J].西安交通大学学报，2014，48（2）：14?19.
XU Shengjun， HAN Jiuqiang， HE Bo， et al. A region Markov random field model with integrated edge feature and image segmentation algorithm [J]. Journal of Xian Jiaotong University， 2014， 48（2）： 14?19.
[7] 余丽玲，徐彬锋，金浩宇，等.基于马尔科夫随机场的乳腺DCE?MRI图像序列配准[J].计算机与现代化，2015（11）：74?78.
YU Liling， XU Binfeng， JIN Haoyu， et al. Images registration of breast DCE?MRI by Markov random field [J]. Computer and modernization， 2015（11）： 74?78.
[8] 林江，戴齐，欧阳婷雪，等.一种边界和马尔可夫随机场相结合的脑MRI医学图像分割方法[J].中国医学物理学杂志，2015，32（5）：717?720.
LIN Jiang， DAI Qi， OUYANG Tingxue， et al. Brain magnetic resonance image segmentation combined boundary and Markov random field [J]. Chinese journal of medical physics， 2015， 32（5）： 717?720.
[9] 郝培博，陈震，江少锋，等.基于Demons算法的2D、3D多模态医学图像非刚性配准研究[J].生物医学工程学杂志，2014（1）：161?165.
HAO Peibo， CHEN Zhen， JIANG Shaofeng， et al. Research on non?regid registration of multi?model medical image based on demons algorithm [J]. Journal of biomedical engineering， 2014（1）： 161?165.
[10] 刘晓慧，杨新锋.基于归一化互信息和金字塔分解优化的图像配准算法研究[J].微型电脑应用，2016，32（11）：19?22.
LIU Xiaohui， YANG Xinfeng. Research on image registration algorithm based on normalized mutual information and pyramid decomposition optimization [J]. Microcomputer applications， 2016， 32 （11）： 19?22.
[11] 秦绪佳，肖佳吉，陈珊，等.局部更新的分层B样条医学图像非刚性配准算法[J].小型微型计算机系统，2016，37（10）：2338?2342.
QIN Xujia， XIAO Jiaji， CHEN Shan， et al. Local updating algorithm of hierarchical B?spline based non?rigid registration for medical images [J]. Journal of Chinese computer systems， 2016， 37（10）： 2338?2342.
摘 要： 在医学图像处理领域中，医学图像配准技术极其重要，其价值体现在临床医学中对图像处理技术的应用。在解决多模态图像配准的相关问题时，基于互信息方法的应用最广泛，但在某些特定的应用中该方法受到的约束仍然较多。针对这一情况，提出一种新的医学图像配准算法，模态变换的引入作为此研究算法的基础，之后新的马尔可夫能量函数则根据两幅通过模态变换后的图像矩阵以及原配准图像得以构建。同时，为了优化能量函数引入了一种改进的梯度下降算法，从而得到配准结果。最后，运用不同的医学图像进行配准实验来验证该算法，通过实验证明该配准算法具有良好的有效性及抗噪性能。
关键词： 图像配准； 模态变换； 马尔可夫随机场； 梯度下降算法； 非刚体配准； 图像矩阵
中图分类号： TN911.73?34； TP311.17 文献标识码： A 文章编号： 1004?373X（2018）01?0057?05
Abstract： In the field of medical image processing， the medical image registration technology is extremely important. Its value is embodied in the application of image processing technology in clinical medicine. For solving the relevant issues of multimodal image registration， the method based on mutual information is most?widely used， but constrained in some specific applications. In view of this situation， a new medical image registration algorithm is proposed. The modal transformation is introduced as the basis of the studied algorithm. After that， the new Markov energy function is constructed according to the two image matrixes after modal transformation and original registration image. In order to optimize the energy function， an improved gradient descent algorithm is introduced to get the registration results. The algorithm was verified with registration experiments of different medical images. The experimental results show that the registration algorithm has perfect effectiveness and high anti?noise performance.
Keywords： image registration； modal transformation； Markov random field； gradient descent algorithm； non?rigid registration； image matrix
0 引 言
图像配准是将图像信息进行汇总并对齐图像的过程，多模态图像配准即为一种将多模态图像进行空间几何关系对齐的技术。近年来，将临床医学和图像处理技术相结合，更加促进了医学图像配准技术的发展。临床上的大多数需求可以被目前已有的算法满足，但是考虑到图像个体之间内容的差异性比较大，而且包含的信息也比较多，所以现今的各种配准算法中都有一定的限制因素存在，一般只能在特定情况下满足配准。考虑到在临床医学方面图像配准具有广泛的应用，而每一种不同的应用均需要不同的形变场，因此关于医学图像配准的问题并没有一种具有普适性的算法解决。针对这一情况，本文提出了一种新的医学图像配准算法。
1 MRF随机场的组成
在图像分析中，马尔可夫随机场的问题可以视为一个建模问题，具体到图像配准中即构造一个有待优化的能量函数问题。基于马尔可夫随机场构造的能量函数（或称为目标函数）不仅和其相邻的像素点有关，而且还需要像素点之间的其他信息，这些信息可以视为一个约束项。一个能量函数可以从以下步骤推导得出。
定义[I（x）]和[J（x）]是二维空间中两幅大小相同的图像，[T（·）]是每一个像素点[x]从固定图像[I（x）]到[J（x）]的一个变换模型，对于两幅给定的图像来说，定义一个观测模型如下所示：
[I（x）=JT（x）+n（x）] （1）
式中：[n（x）]是一个高斯白噪声模型，其对应的噪声模型数学期望为[x～（0，σ2）]，[σ]是噪声模型的标准差。
将两幅图像和变换模型视为随机变量，[s]和[t]表示像素位置，则可以通过图[?=（v，ε）]表示马尔可夫随机场模型中离散标记的问题。定义[?=（1，2，…，L）]为图像的一个标号集合，则所有的标号集合可以定义为：
[X=xss∈V] （2）
式中：[xs]是[?]中的一个元素。给定一个包含平移、翻转、旋转和在[x，y]方向上的尺度变换的向量[θ，]可以得到标号集合的马尔可夫随机场模型的能量函数：
[Emrf（Xθ）=Edata+Esmooth] （3）
1.1 数据项
式（3）反映了两个相对应的像素点及其邻域像素点之间的关系。通常来说，在整个定义域[Ω]中，构造一个解决配准问题的模型是定义一个数据项。这个数据项能够表示在固定图像和浮动图像之间的距离函数。根据式（1）和式（2），图像[I（x）]的噪声模型的期望是[JT（x）]，方差是[σ2，]因此马尔可夫随机场模型的数据项可以表示为：
[Edata=s∈Vθs（xs）=12σ2s∈VI（xs）-JT（x）2] （4）
参数[σ]是作为测度项的调节参数引入的，它的具体值取决于人为选定的噪声模型以及实际的应用。尽管只采用数据项进行优化的方法在某些特定应用中是可行的，但是这并非一种十分通用的方法。这一方法的局限性体现在：某些应用中方程的解不存在；某些情况下对该问题的解不惟一；某些问题中，方程的解对原始数据的依赖不具有连续性。
1.2 平滑项
上述提到的所有限制都可以归结成一个不适定问题（ill?posed problem），这一问题最直接的解决方法是采用一个平滑项做约束。对于低层级的马尔可夫随机场问题，这一项可以视为一个归一化项（regularization）。它假定了在相同的邻域系统中相邻像素点的物理特性是不会突然发生变化的。因此，对于一个邻域系统[Ni]而言，其平滑项为：
[Esmooth=λs，t∈Niθstxs，xt=λs∈Vt∈Niθstxs，xt] （5）
式中：[λ]表示一个权重值，其值通常根据不同的配准形式以及配准所采用的不同函数[θst]来选取不同的常数值。函数[θst]的选择有很多种，比如分段函数、多项式函数、对数函数等。鉴于这些函数在某些情况下会加大计算的复杂度，本文采用一种线性函数来简化优化过程。采用一种间断自适应（Discontinuity Adaptive， DA）模型来构造函数的平滑项，而且函数[θ（x）]需要满足以下方程式：
[limx→∞θ（x）=limx→∞2xf（x）=C] （6）
式中：[C]是一个正整数；函数[f（x）]是一个在定义域中的偶函数。一般来说，马尔可夫随机场模型通常是指二阶的随机场模型。对于一个给定的概率模型，马尔可夫随机场问题可以看成是一个解决最大后验概率的问题。在上述讨论条件不变的情况下将固定图像、浮动图像和变换模型视为随机变量，则根据三者之间的统计相关性可得：
[pI（xs）T，J（xt）=pIJT（xt）] （7）
但是对于一个随机变量[X，]上述的分布往往难以求得，此时，可以考虑运用Hammersley?Clifford 定理将马尔可夫随机场模型转换为吉布斯分布：
[p（Xθ）∝e-E（Xθ）] （8）
根据式（8）可得，将马尔可夫的最大后验概率的相关约束问题换一种思路进行考虑，即用最小化能量函数[minXE（Xθ）]的相关约束问题来代替则更容易解决。
2 需求分析和系统设计
为了更大限度地提高医学图像配准的精度，本文在对目标函数进行优化之前引入了一个基于直方图计算的模态变换方法。
假设有两幅不同模态的图像[I（x）]和[J（x）]已经归一化，同时[I（x）]，[J（x）]∈[0，1]，[x]为像素点位置，通过对每一个像素点进行一个遍历，每遍历一个点就更新一次联合直方图，直到最后一个点遍历完，得到图像[I]和[J]的联合直方图[H（I，J）]。
[HI（x）N，J（x）N=HI（x）N，J（x）N+1] （9）
式中：[N]是联合直方图的样条数；[x]是当前像素点的位置；[c]表示小于[c]的最大正整数值；联合直方图[H（i，j）]是一个二维矩阵；灰度值[i]在图[I]以及灰度值[j]在图像[J]中相关点具体的数目用[H（i，j）]来表示。
运用式（10）和式（11）能够将图像[I]和图像[J]的灰度值表示形式的图像矩阵转换得大致相似。通过迭代过程，搜索图像[I]在图像[J]中重叠次数最多的一个像素点。在遍历了所有像素点之后，得到一个新的图像矩阵，也即经过模态变换后的图像矩阵。同理，可以由图像[J]得到[JT。]
[IT（x）=argmaxjHI（x）N，jN] （10）
[JT（x）=argmaxiHiNN，J（x）N] （11）
图1给出了利用MRI?T1序列脑部图像和CT图像进行模态变换的结果图。图1c）和图1d）是进行模态变换前后的联合直方图，从图中可以明显看到在进行模态变换前，联合直方图上的点分布比较松散，变换之后的点分布紧凑，且向图中的斜对角线靠拢，这表明模态变换的方法减少了图像簇的分散性，在某种程度上可以提高图像的配准精度。
3 基于MRF的模态变换多模态图像配准
3.1 MRF能量函数
考虑到图像[I]以及图像[J，]在图像的定义域[Ω=1，m×1，n]已经给定的情况下，将得到的4个图像矩阵结合可以得到如下数据项：
[Edata=12σ21s∈VIT（xs）-JT（x）2+12σ22s∈VI（xs）-JTT（x）2] （12）
式中：[σ1]和[σ2]是高斯白噪声模型的标准差。在实验中，为了简化计算过程，取[σ1=σ2=σ]。通常来说，数据项来自于对距离函数的定义，它通常处理的是大多数的线性问题。为了提高配准精度以及配准过程的稳定性，本文引入了另外一项表达式表示数据项，更多的细节信息的获得來自于之前转换后的图像矩阵[IT]和[JT。]
能量函数只和局部形变有关，但是在优化过程中仍然有可能会产生局部极值的情况，因此本文选择了一个基本函数作为平滑项的势能函数：
[θ（x）=ln（1+x2）] （13）
[f（x）=11+x2] （14）
式中：函数[θ（x）]和[f（x）]均满足表达式（6）的条件。由于在形变过程中，平滑项更多的是和图像像素的邻域系统有关，同时结合表达式（5）、（12）和（13），得到平滑项的相关计算表达式如下：
[Esmooth=λs∈Vt∈Niln1+xs-xt2] （15）
平滑项的引入能够处理非线性问题，同时平滑形变场，从而可以在优化过程中减少极值对配准结果的影响。
在得到数据项和平滑项之后，根据马尔可夫随机场模型理论可以得到一个新的能量函数模型：
[Emrf=Edata+Esmooth=12σ21s∈VIT（xs）-JT（x）2+ 12σ22s∈VI（xs）-JTT（x）2+λs∈Vt∈Niln1+xs-xt2] （16）
值得注意的一点是，当图像完全配准之后，马尔可夫随机场的能量函数值[Emrf]将会最小。通常来说，一个马尔可夫随机场模型指的是一个二阶模型，在计算过程中需要求像素点相邻的八个像素点（一阶系统中，斜对角线上的像素点不考虑在内）。由于二阶邻域系统计算了更多的图像信息，从而可以得到更精确的配准结果，本文中的实验均采用二阶邻域系统。
3.2 优化算法
对于医学图像配准来说，优化方法主要分为两大类：参数化和非参数化配准。因为非参数化配准广泛应用于非刚体医学图像配准中，本文主要讨论这一种优化方法。优化方法的通用模型如下：
[xi+1=xi+ai?di，i=0，1，2，…] （17）
式中：[xi]和[xi+1]分别为第[i]次和第[i+1]次迭代过程中的优化值；[di]表示当前优化过程的搜索方向，搜索方向既可以为正值也可以为负值。沿着搜索方向，[ai]是一个在每次迭代过程中用来控制步长的标量尺度系数。
一个最基本的优化方法是梯度下降方法，它是众多无约束最优化方法中的一种。这种方法迭代简单，而且相对应的配准结果也较为准确。但是这种算法在优化时并不是逐次向最优项靠近，而且由于迭代步长较小，所以迭代时间相对较长。因此，通过用图形的子集计算互信息的差分方法来提高多模态图像配准的计算时间。这一方法用预定义的一个衰减函数取代了常数值作为每次迭代的增益。其中衰减函数的表达式如下：
[ai=ai+Aα] （18）
式中：[a>0，A≥1，0<α≤1。]这一函数能够更好地选择步长值，同时也减少了优化过程中的迭代次数。在此基础上，用自适应随机梯度算法（Adaptive Stochastic Gradient Descent，ASGD）对该算法进行改进。这一方法通过计算[di+1]和[di]内积值来调整步长[ai]的大小：
[xi+1=xi-ai?di，i=0，1，2，…] （19）
[ai+1=ai?f-dTi?di-1+] （20）
式中：[[x]+]表示[x]和0进行比较后取较大的数，即[max（x，0）；]非线性sigmoid函数用函数[f]进行表示，该函数的一般表达式为[f（x）=11+e-xw]。在迭代过程中，若目标函数的当前值小于最后一次迭代值，则保持当前步长不变，否则当前步长将被式（20）的步长值取代。
通过将上述两种优化方法进行结合的方法改进优化算法。在同一个配准实验中，这一方法提供了一个更加灵活的选择步长的方式，同时使得在优化过程中减少了迭代次数，且配准输出的结果和文献中的算法具有可比性。算法中的步长值[ai+1]定义为：
[ai+1=ai?f-dTi?di-1+，Δ≥εα?ai，Δ<ε] （21）
式中：[α]是一个在定义域为[0.9，1]的常数；[Δ]是[ai]和[ai+1]差值的绝对值[ai-ai+1]；[ε]是一个根据经验而选取的值，通常它是一个小于1的常数。将梯度[di]由[ΔEmrf]代替，即得到最终的优化算法：
[xi+1=xi-ai?ΔEmrf（xi）] （22）
由式（22）可以看到，对于能量函数以及能量函数的导数，通常需要逐点进行求解，这样对于较大尺寸的图像进行配准时，会相当耗时。因此，文中所做的实验均采用一种多分辨率（multiresolution）的配准方法。这一方法的引入可以避免配准过程中遇到的局部极值，从而提高配准精度。将这一方法应用到研究中，首先对一个16×16的子网格像素点进行配准，然后将得到的形变场和原始图像调整到32×32的子网格再进行配准，以此类推，直到配准过程中这一网格达到待配准图像的原始大小，配准结束。图2中，三条不同的曲线分别表示在每一个层级优化过程的迭代情况。
为了验证本文优化方法的有效性，本文将上述方法和前人所用的研究方法进行了对比。实验中除了优化方法以外，所有的參数都设定为相同值。在式（18）和式（21）中，参数[A=10，a=1，α=0.6，ε=0.1。]同时，多分辨率方法设定层次为三层。实验用脑部图像参见图1。迭代过程中两种优化方法在迭代过程中的曲线如图2所示。从图2a）中可以看出，Myronenko方法在前20～40次迭代时陷入了局部极值，且这一情况在优化的最后一个层级仍然存在。而本文提出的优化算法不仅在迭代次数上少于Myronenko方法（本文算法的迭代次数少于100次），而且在整个迭代过程中曲线较为平滑，较大程度上解决了图2a）出现的局部极值问题。
4 算法验证
为了验证本文所提算法的可行性以及配准效率，本研究做了一系列图像配准实验。当配准过程开始时，需要应用一个全局性的仿射变换以增加配准的准确性。初始的仿射变换参数[x=][0，0，0，100，100，0，0]分别代表在[x，y]方向上的平移、旋转、尺度变换以及裁剪（shearing）。之后，一个三层的多分辨率方法同时应用于固定图像和浮动图像，所有图像的灰度值均归一化为[0，1]的区间。在每一个层级的迭代过程中，都应用模态变化的方法。之后，优化方法对已经构建好的能量函数进行优化，能量函数的参数设定为[σ1=σ2=σ=1]。所有最大迭代次数都设定为200次。
在第一个刚体实验中，选取了MRI?T1序列脑部图像和MRI?PD序列脑部图像分别作为固定图像和浮动图像，可参见图1a）和图1b）。浮动图像是在固定图像的基础上人为旋转了10°，同时在[x]和[y]方向上分别有13 mm和17 mm的位移。图像大小均为221×257。图3给出了两幅实验图像的配准后图像以及配准前后的棋盘式叠加图。通过棋盘图可以看到图像边缘有了很好的配准效果。
实验2中，验证了3种不同的非刚体形变情况下，本文算法的配准效果。如图4所示为本实验的配准用图及配准结果。图中的每一列分别为固定图像、浮动图像、配准后图像以及配准前后的棋盘格叠加图。从配准前后叠加图的效果来看，本文提出的算法输出了准确的配准结果，图中从上到下的图像模态分别为MRI?T2/CT、MRIT2/MRIT1和CT/MRI?T2。
5 结 语
本文在结合MRF随机场和模态变换的基础上，构造了一个新的待优化能量函数。此能量函数由数据项和平滑项构成。数据项包含了4个图像矩阵，分别为固定图像和浮动图像，以及由这两幅图像进行模态变换后得到的2个图像矩阵。增加两幅图像矩阵可以增加图像的细节，通过获得更多的图像信息从而提高图像配准的精度，平滑项选取了自然对数函数。
在优化方法的选取上，本文改进了梯度下降算法，将两种优化方法进行结合，使优化过程中的迭代步长选择更加灵活。通过实验验证，此方法能够有效地减少图像迭代次数，大大降低了优化过程陷入局部极值的情况。为了验证本文算法的配准精度以及抗噪性能，本文用MRI不同序列脑部图像和CT图像做了一系列的配准实验。实验结果表明，本文提出的配准方法配准效果更为准确，同时抗噪性能更好。
参考文献
[1] 张冉，王雷，夏威，等.2D/3D图像配准中的相似性测度和优化算法[J].激光与红外，2014（1）：98?102.
ZHANG Ran， WANG Lei， XIA Wei， et al. Comparison of similarity measurement and optimization [J]. Laser & infrared， 2014（1）： 98?102.
[2] 王雷，黄晨雪，王高波.改进的MRF彩色图像分割算法[J].计算机工程与应用，2016，52（13）：191?194.
WANG Lei， HUANG Chenxue， WANG Gaobo. Improved MRF color image segmentation algorithm [J]. Computer engineering and applications， 2016， 52（13）： 191?194.
[3] 张静亚，王加俊.一种改进的非刚性医学图像配准算法[J].计算机应用研究，2015，32（4）：1261?1264.
ZHANG Jingya， WANG Jiajun. Improved non?rigid medical image registration algorithm [J]. Application research of computers， 2015， 32（4）： 1261?1264.
[4] 金斌，周伟，丛瑜，等.一种基于局部不变特征的SAR图像配准新算法[J].哈尔滨工业大学学报，2014，46（11）：112?118.
JIN Bin， ZHOU Wei， CONG Yu， et al. A novel and efficient algorithm using local invariant feature for SAR image registration [J]. Journal of Harbin Institute of Technology， 2014， 46（11）： 112?118.
[5] 赵雪梅，李玉，赵泉华.结合高斯回归模型和隐马尔可夫随机场的模糊聚类图像分割[J].电子与信息学报，2014，36（11）：2730?2736.
ZHAO Xuemei， LI Yu， ZHAO Quanhua. Image segmentation by fuzzy clustering algorithm combining hidden Markov random field and Gaussian regression model [J]. Journal of electronics & information technology， 2014， 36（11）： 2730?2736.
[6] 徐胜军，韩九强，何波，等.融合边缘特征的马尔可夫随机场模型及分割算法[J].西安交通大学学报，2014，48（2）：14?19.
XU Shengjun， HAN Jiuqiang， HE Bo， et al. A region Markov random field model with integrated edge feature and image segmentation algorithm [J]. Journal of Xian Jiaotong University， 2014， 48（2）： 14?19.
[7] 余丽玲，徐彬锋，金浩宇，等.基于马尔科夫随机场的乳腺DCE?MRI图像序列配准[J].计算机与现代化，2015（11）：74?78.
YU Liling， XU Binfeng， JIN Haoyu， et al. Images registration of breast DCE?MRI by Markov random field [J]. Computer and modernization， 2015（11）： 74?78.
[8] 林江，戴齐，欧阳婷雪，等.一种边界和马尔可夫随机场相结合的脑MRI医学图像分割方法[J].中国医学物理学杂志，2015，32（5）：717?720.
LIN Jiang， DAI Qi， OUYANG Tingxue， et al. Brain magnetic resonance image segmentation combined boundary and Markov random field [J]. Chinese journal of medical physics， 2015， 32（5）： 717?720.
[9] 郝培博，陈震，江少锋，等.基于Demons算法的2D、3D多模态医学图像非刚性配准研究[J].生物医学工程学杂志，2014（1）：161?165.
HAO Peibo， CHEN Zhen， JIANG Shaofeng， et al. Research on non?regid registration of multi?model medical image based on demons algorithm [J]. Journal of biomedical engineering， 2014（1）： 161?165.
[10] 刘晓慧，杨新锋.基于归一化互信息和金字塔分解优化的图像配准算法研究[J].微型电脑应用，2016，32（11）：19?22.
LIU Xiaohui， YANG Xinfeng. Research on image registration algorithm based on normalized mutual information and pyramid decomposition optimization [J]. Microcomputer applications， 2016， 32 （11）： 19?22.
[11] 秦绪佳，肖佳吉，陈珊，等.局部更新的分层B样条医学图像非刚性配准算法[J].小型微型计算机系统，2016，37（10）：2338?2342.
QIN Xujia， XIAO Jiaji， CHEN Shan， et al. Local updating algorithm of hierarchical B?spline based non?rigid registration for medical images [J]. Journal of Chinese computer systems， 2016， 37（10）： 2338?2342.
