CN109741378A - 基于mrf模型的多模态医学图像配准方法、装置、平台及介质 - Google Patents

Abstract

本发明公开了一种基于MRF模型的多模态医学图像配准方法，具体步骤包括：输入待配准的两幅图像，分别记为固定图和浮动图；构建马尔可夫随机场模型的能量函数；将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；对最终得到的配准结果进行显示。本发明引入了模态变换，并根据两幅通过模态变换后的图像矩阵以及原配准图像来构建新的马尔可夫能量函数，同时还通过引入一种改进的梯度下降算法来优化能量函数，从而得到配准结果。最后通过不同的医学图像进行配准实验来验证本发明，证明了本发明具有良好的有效性及抗噪性能。

Description

基于MRF模型的多模态医学图像配准方法、装置、平台及介质

技术领域

本发明涉及医学图像处理领域，尤其涉及一种基于MRF模型的多模态医学图像配准方法、装置、平台及介质。

背景技术

图像配准是将图像信息进行汇总并对齐图像的过程，而多模态图像配准是一种将多模态图像进行空间几何关系对齐的技术。目前，临床上的大多数需求可以被目前已有的算法满足，但是考虑到图像个体之间内容的差异性比较大，而且包含的信息也比较多，所以现今的各种配准算法中都有一定的限制因素存在，一般只能在特定情况下满足配准。考虑到在临床医学方面图像配准具有广泛的应用，而每一种不同的应用均需要不同的形变场，因此关于医学图像配准的问题并没有一种具有普适性的算法解决。

目前常见的图像配准技术包括基于改进互信息和相关率的图像配准方法、基于特征点和图像分割的配准方法以及基于熵图和拉普拉斯图的图像配准方法等。其中，基于互信息的图像配准方法相比于其他图像配准方法，计算复杂度高，用时长。而基于特征点和图像分割的配准技术首先要确定标记点或分割图像，但由于超声图像质量不高，使得自动选取标记点或者图像分割的误差都比较大，人工选取标记点或人工图像分割的实时性不高。对于基于特征描述的图像配准方法，其基于图像块的运算会丢失部分局部信息，影响配准精度，且该类方法适用于较小变形，随着变形程度的增加，配准误差增大。

综上所述的图像配准技术在解决多模态图像配准的相关问题时，基于互信息方法的应用最广泛，但在某些特定的应用中该方法受到的约束仍然较多。

发明内容

本发明的第一目的在于克服现有技术的不足，提供一种基于MRF模型的多模态医学图像配准方法，本发明以模态变换为基础，根据两幅经过模态变换后的图像矩阵以及原配准图像构建新的马尔可夫能量函数，同时采用改进的梯度下降算法来优化能量函数，从而得到配准结果。本发明具有良好的有效性以及抗噪性。

本发明的第二目的在于提供一种基于MRF模型的多模态医学图像配准装置。

本发明的第三目的在于提供一种平台。

本发明的第四目的在于提供一种存储介质。

本发明的第一目的能够通过以下技术方案实现：

一种基于MRF模型的多模态医学图像配准方法，具体步骤包括：

输入待配准的两幅图像，分别记为固定图和浮动图；

构建马尔可夫随机场模型的能量函数；

将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；

对最终得到的配准结果进行显示。

具体地，所述构建马尔可夫随机场模型的能量函数，包括：

S201、获取输入的待配准图像；

S202、对浮动图进行处理，并将固定图与处理好的浮动图进行归一化；

S203、采用权重图的形式对固定图和浮动图的对应的特殊位置进行标记；

S204、计算固定图和浮动图每个像素之间的距离值，将得到的相似度函数作为数据项；

S205、生成能量函数的平滑项；

S206、根据数据项和平滑项得到马尔可夫随机场模型的能量函数；

S207、将马尔可夫随机场模型的能量函数进行转换；

S208、如果遍历完整图像的所有像素点，则结束；否则转到S204；

S209、根据固定图和处理好的浮动图的归一化结果更新一次联合直方图；

S210、对固定图的一个像素进行迭代，搜索在固定图中重叠次数最多的像素点；

S211、对浮动图的一个像素进行迭代，搜索在浮动图中重叠次数最多的像素点；

S212、判断是否遍历完整图像的所有像素点；

S213、获得模态变换后的固定图与浮动图。

具体地，将两幅图像和变换模型视为随机变量，则可以通过图表示马尔可夫随机场模型中离散标记的问题，定义为图像的一个标号集合，则所有的标号集合可以定义为：

X＝{x_s|s∈V}

其中，s表示标记的像素位置，t表示未做标记的像素位置，V表示s的集合，θ表示t的集合，L表示标记的顺序号，x_s为中的一个元素。

进一步地，步骤S202中对浮动图进行处理包括：平移、翻转、旋转以及在x-y方向上的尺度变换。

进一步地，步骤S204中根据构建的距离函数来表示相似度，表示方式为：

E_data＝∑_s∈Vθ_s(x_s)＝1/2σ²∑_s∈V[I(x_s)-J(T(x))]²

给定图像的定义域，将得到的四个图像矩阵结合得到数据项，表示为：

E_data＝1/2σ₁ ²∑_s∈V[I(x_s)-J(T(x))]²+1/2σ₂ ²∑_s∈V[I(x_s)-J_t(T(x))]²

其中，σ₁和σ₂为高斯白噪声模型的标准差，在实验中为了简化计算过程，取σ₁＝σ₂＝σ。

进一步地，在步骤S205中，采用一种间断自适应模型构造能量函数的平滑项，用于作为约束来解决不适定问题，表示为：

E_smooth＝λ∑_s,t∈Niθ_st(x_s,x_t)＝λ∑_s∈V∑_t∈θθ_st(x_s,x_t)

Ni表示邻域系统。

能量函数只和局部形变有关，但是在优化过程中仍然有可能会产生局部极值的情况，因此本发明选择了一个基本函数作为平滑项的势能函数：

Θ(x)＝ln(1+x²)

F(x)＝1/(1+x²)

平滑项的相关计算表达式表示为：

E_smooth＝λ∑_s,t∈Niln[1+(x_s-x_t)²]

平滑项的引入能够处理非线性问题，同时平滑形变场，从而可以在优化过程中减少极值对配准结果的影响。

进一步地，所述步骤S206中得到的能量函数表示为：

E_mrf(X|θ)＝E_data+E_smooth

在根据步骤S204和S205得到的数据项和平滑项之后，根据马尔可夫随机场模型理论可以得到一个具体的能量函数模型，表示为：

E_mrf＝E_data+E_smooth＝1/2σ₁ ²∑_s∈V[I(x_s)-J(T(x))]²

+1/2σ₂ ²∑_s∈V[I(x_s)-J_t(T(x))]²+λ∑_s,t∈Niln[1+(x_s-x_t)²]

当图像完全配准之后，马尔可夫随机场的能量函数值将会最小。

进一步地，在步骤S207中运用Hammersley Clifford定理将模型转换为吉布斯分布，表示为：

进而转为最小化能量函数的相关约束问题。

进一步地，在步骤S209中，假设有两幅不同模态的图像I(x)和J(x)已经归一化，同时I(x),J(x)∈[0，1]，x为像素点位置，通过对每一个像素点进行一个遍历，每遍历一个点就更新一次联合直方图，直到最后一个点遍历完，得到图像I和J的联合直方图H(I,J)；联合直方图更新的公式为：

其中，N表示联合直方图的样条数；x表示当前像素点的位置；表示小于c的最大正整数值；联合直方图H(i,j)是一个二维矩阵；灰度值i在图I(x)以及灰度值j在图像J(x)中相关点具体的数目用H(i,j)来表示。

通过迭代过程，搜索图像I(x)在图像J(x)中重叠次数最多的一个像素点。在遍历了所有像素点之后，得到一个新的图像矩阵，也即经过模态变换后的图像矩阵。同理，可以由图像J(x)得到J_T(x)。具体转换公式如下所示：

具体地，所述将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数，包括：

S301、初始化sigmoid函数以及优化模块的常数；

S302、初始化步长a₀，根据公式计算出a₁；

S303、根据最终的能量函数，计算固定图和浮动图每个像素及搜索方向；

S304、对固定图和浮动图的像素点分别进行迭代更新；

S305、如果x_i+1＜x_i成立，转到S310；否则转到S306；

S306、如果Δ＝|a_i-a_i-1|≥ε转到S307；否则转到S308；

S307、更新a_i+1＝[a_i·f(-d_i ^T·d_i-1)]⁺；

S308、更新a_i+1＝α·a_i；

S309、令i＝i+1，重新转到S304；

S310、如果遍历完图像所有像素，转到S311；否则转到S304；

S311、记录配准得到的形变场结果、原始图像以及模态变换后的图像，流程结束。

进一步地，在步骤S301中，采用预定义的一个衰减函数取代常数值作为每次迭代的增益，衰减函数的表达式为：

a_i＝a/(i+A)^α

其中，a＞0，A≥1，0＜α≤1，A表示衰减常数。这一函数能够更好地选择步长值，同时减少了优化过程中的迭代次数。在此基础上，用自适应随机梯度算法对该算法进行改进，具体为：

x_i+1＝x_i-a_i·d_i,i＝0,1,2,...

a_i+1＝[a_i·f(-d^T·d_i-1)]⁺

其中，i表示迭代次数，d_i表示第i次迭代过程的搜索方向，搜索方向既可以为正值也可以为负值，[x]⁺表示x和0进行比较后取较大的数，即max(x,0)；非线性sigmoid函数用函数f进行表示，其一般表达式为f(x)＝1/(1+e-x/w)；其中，d^T表示方向矩阵的转置、e为常数、w表示输入向量、x表示图像的像素。在迭代过程中，若目标函数的当前小于最后一次迭代值，则保持当前步长不变，否则当前步长将被上式的步长值取代。

将上述两种优化方法进行结合来改进优化算法，改进后算法中步长值定义为：

其中，α表示一个定义域为[0.9，1.0]的常数，Δ表示a_i和a_i+1差值的绝对值，即|a_i-a_i+1|；ε是一个根据经验选取的值，通常为一个小于1的常数。将d_i采用ΔE_mrf进行替代，得到最终优化算法为：

x_i+1＝x_i-a_i·ΔE_mrf(x_i)

具体地，所述对配准结果进行显示的步骤中，采用一种多分辨率的配准方法显示最终的配准结果，配准过程为：先对一个16×16的子网格像素点进行配准，将得到的形变场和原始图像调整到32×32的子网格再进行配准，以此类推，当配准的网格达到待配准图像的原始大小时配准结束，对配准结果进行显示。

本发明的第二目的能够通过以下技术方案实现：

一种基于MRF模型的多模态医学图像配准装置，包括：

图像生成模块，用于输入待配准的两幅图像，分别记为固定图和浮动图；

能量函数构建模块，用于构建马尔可夫随机场模型的能量函数；

最优化计算模块，用于将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；

配准结果显示模块，用于对最终得到的配准结果进行显示。

进一步地，在所述能量函数构建模块中，包括：

图像转化单元，用于将浮动图进行任意的平移、翻转、旋转以及在x-y方向上的尺度变换；

图像归一化单元，用于将固定图以及浮动图进行归一化处理；

离散标记图生成单元，用于记录固定图和浮动图对应特殊位置的标记，并转化为权重图的形式；

数据项生成单元，用于通过计算固定图和浮动图每个像素之间的距离值，得到相似度函数作为数据项；

平滑项生成单元，用于生产能量函数的平滑项，作为函数模型的约束项；

能量函数生成单元，用于根据平滑项和数据项生产对应的MRF能量函数；

能量函数转化单元，用于将模型转化为最小化能量函数的相关约束；

模态变换单元，用于对固定图和浮动图进行一次基于直方图计算的模态变换；

进一步地，在所述最优化计算模块中，包括：

初始化单元，用于初始化sigmoid函数、优化模块的常数以及初始化步长a₀，并根据步长公式计算出a₁；

像素计算单元，用于根据最终的能量函数，计算固定图和浮动图每个像素及搜索方向；

更新单元，用于对固定图和浮动图的像素点以及步长分别进行迭代更新；

判断单元，用于对最优化计算过程中参数进行迭代更新的条件进行判断；

记录单元，用于记录配准得到的形变场结果、原始图像以及模态变换后的图像。

本发明的第三目的能够通过以下技术方案实现：

一种平台，包处理器以及用于存储处理器可执行程序的存储器，所述处理器执行存储器存储的程序时，实现上述的图像配准方法。

本发明的第四目的能够通过以下技术方案实现：

一种存储介质，存储有程序，所述程序被处理器执行时，实现上述的图像配准方法。

本发明相较于现有技术，具有以下的有益效果：

1、本发明中所使用的形变场更具普遍性。

2、本发明能够有限减少图像迭代次数，降低优化过程陷入局部极值的情况。

3、本发明得到的配准效果更为准确，同时具有更好的抗噪性能。

附图说明

图1为本发明实施例中的基于MRF模型的多模态医学图像配准方法的流程图。

图2为本发明实施例中的构建马尔可夫随机场模型的能量函数的流程图。

图3为本发明实施例中的采用迭代方式，将模型输入到改进的优化模型进行最优化计算，得到最终变换参数的流程图。

图4为本发明实施例中的基于MRF模型的多模态医学图像配准装置的结构示意图。

图5为本发明实施例中运用配准方法后图像示意图和运用配准方法前后图像的叠加图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1：

如图1所示，本实施例提供了一种基于MRF模型的多模态医学图像配准方法，该方法具体步骤包括：

输入待配准的两幅图像，分别记为固定图和浮动图；

构建马尔可夫随机场模型的能量函数；

将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；

对最终得到的配准结果进行显示。

如图2所示为构建马尔可夫随机场模型的能量函数步骤的具体流程图，包括：

S201、获取输入的待配准图像；

S202、对浮动图进行处理，并将固定图与处理好的浮动图进行归一化；

S203、采用权重图的形式对固定图和浮动图的对应的特殊位置进行标记；

S204、计算固定图和浮动图每个像素之间的距离值，将得到的相似度函数作为数据项；

S205、生成能量函数的平滑项；

S206、根据数据项和平滑项得到马尔可夫随机场模型的能量函数；

S207、将马尔可夫随机场模型的能量函数进行转换；

S208、如果遍历完整图像的所有像素点，则结束；否则转到S204；

S209、根据固定图和处理好的浮动图的归一化结果更新一次联合直方图；

S210、对固定图的一个像素进行迭代，搜索在固定图中重叠次数最多的像素点；

S211、对浮动图的一个像素进行迭代，搜索在浮动图中重叠次数最多的像素点；

S212、判断是否遍历完整图像的所有像素点；

S213、获得模态变换后的固定图与浮动图。

进一步地，步骤S202中对浮动图进行处理包括：平移、翻转、旋转以及在xy方向上的尺度变换。

进一步地，步骤S204中根据构建的距离函数来表示相似度，表示方式为：

E_data＝∑_s∈Vθ_s(x_s)＝1/2σ²∑_s∈V[I(x_s)-J(T(x))]²

给定图像的定义域，将得到的四个图像矩阵结合得到数据项，表示为：

E_data＝1/2σ₁ ²∑_s∈V[I(x_s)-J(T(x))]²+1/2σ₂ ²∑_s∈V[I(x_s)-J_t(T(x))]²

其中，σ₁和σ₂为高斯白噪声模型的标准差，在实验中为了简化计算过程，取σ₁＝σ₂＝σ。

进一步地，在步骤S205中，采用一种间断自适应模型构造能量函数的平滑项，用于作为约束来解决不适定问题，表示为：

E_smooth＝λ∑_s,t∈Niθ_st(x_s,x_t)＝λ∑_s∈V∑_t∈θθ_st(x_s,x_t)

能量函数只和局部形变有关，但是在优化过程中仍然有可能会产生局部极值的情况，因此本发明选择了一个基本函数作为平滑项的势能函数：

Θ(x)＝ln(1+x²)

F(x)＝1/(1+x²)

平滑项的相关计算表达式表示为：

E_smooth＝λ∑_s,t∈Niln[1+(x_s-x_t)²]

平滑项的引入能够处理非线性问题，同时平滑形变场，从而可以在优化过程中减少极值对配准结果的影响。

进一步地，所述步骤S206中得到的能量函数表示为：

E_mrf(X|θ)＝E_data+E_smooth

在根据步骤S204和S205得到的数据项和平滑项之后，根据马尔可夫随机场模型理论可以得到一个具体的能量函数模型，表示为：

E_mrf＝E_data+E_smooth＝1/2σ₁ ²∑_s∈V[I(x_s)-J(T(x))]²

+1/2σ₂ ²∑_s∈V[I(x_s)-J_t(T(x))]²+λ∑_s,t∈Niln[1+(x_s-x_t)²]

当图像完全配准之后，马尔可夫随机场的能量函数值将会最小。

进一步地，在步骤S207中运用Hammersley Clifford定理将模型转换为吉布斯分布，表示为：

进而转为最小化能量函数的相关约束问题。

如图3所示为将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数步骤的具体流程图，包括：

S301、初始化sigmoid函数以及优化模块的常数；

S302、初始化步长a₀，根据公式计算出a₁；

S303、根据最终的能量函数，计算固定图和浮动图每个像素及搜索方向；

S304、对固定图和浮动图的像素点分别进行迭代更新；

S305、如果x_i+1＜x_i成立，转到S310；否则转到S306；

S306、如果Δ＝|a_i-a_i-1|≥ε转到S307；否则转到S308；

S307、更新a_i+1＝[a_i·f(-d_i ^T·d_i-1)]⁺；

S308、更新a_i+1＝α·a_i；

S309、令i＝i+1，重新转到S304；

S310、如果遍历完图像所有像素，转到S311；否则转到S304；

S311、记录配准得到的形变场结果、原始图像以及模态变换后的图像，流程结束。

进一步地，在将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数的步骤中，可以采用预定义的一个衰减函数取代常数值作为每次迭代的增益，衰减函数的表达式为：

a_i＝a/(i+A)^α

其中，a＞0，A≥1，0＜α≤1。这一函数能够更好地选择步长值，同时减少了优化过程中的迭代次数。在此基础上，用自适应随机梯度算法对该算法进行改进，具体为：

x_i+1＝xi-a_i·d_i,i＝0,1,2,…

a_i+1＝[a_i·f(-d^T·d_i-1)]⁺

其中，[x]⁺表示x和0进行比较后取较大的数，即max(x,0)；非线性sigmoid函数用函数f进行表示，其一般表达式为f(x)＝1/(1+e-x/w)。在迭代过程中，若目标函数的当前小于最后一次迭代值，则保持当前步长不变，否则当前步长将被上式的步长值取代。

将上述两种优化方法进行结合来改进优化算法，改进后算法中步长值定义为：

其中，α表示一个定义域为[0.9，1.0]的常数，Δ表示a_i和a_i+1差值的绝对值，即|a_i-a_i+1|；ε是一个根据经验热选取的值，通常为一个小于1的的常数。将d_i采用ΔE_mrf进行替代，得到最终优化算法为：

x_i+1＝x_i-a_i·ΔE_mrf(x_i)

具体地，所述对配准结果进行显示的步骤中，采用一种多分辨率的配准方法显示最终的配准结果，配准过程为：先对一个16×16的子网格像素点进行配准，将得到的形变场和原始图像调整到32×32的子网格再进行配准，以此类推，当配准的网格达到待配准图像的原始大小时配准结束，对配准结果进行显示。

实施例2：

如图4所示，本实施例一种基于MRF模型的多模态医学图像配准装置，包括：

图像生成模块，用于输入待配准的两幅图像，分别记为固定图和浮动图；

能量函数构建模块，用于构建马尔可夫随机场模型的能量函数；

最优化计算模块，用于将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；

配准结果显示模块，用于对最终得到的配准结果进行显示。

进一步地，在所述能量函数构建模块中，包括：

图像转化单元，用于将浮动图进行任意的平移、翻转、旋转以及在x,y方向上的尺度变换；

图像归一化单元，用于将固定图以及浮动图进行归一化处理；

离散标记图生成单元，用于记录固定图和浮动图对应特殊位置的标记，并转化为权重图的形式；

数据项生成单元，用于通过计算固定图和浮动图每个像素之间的距离值，得到相似度函数作为数据项；

平滑项生成单元，用于生产能量函数的平滑项，作为函数模型的约束项；

能量函数生成单元，用于根据平滑项和数据项生产对应的MRF能量函数；

能量函数转化单元，用于将模型转化为最小化能量函数的相关约束；

模态变换单元，用于对固定图和浮动图进行一次基于直方图计算的模态变换；

进一步地，在所述最优化计算模块中，包括：

初始化单元，用于初始化sigmoid函数、优化模块的常数以及初始化步长a₀，并根据步长公式计算出a₁；

像素计算单元，用于根据最终的能量函数，计算固定图和浮动图每个像素及搜索方向；

更新单元，用于对固定图和浮动图的像素点以及步长分别进行迭代更新；

判断单元，用于对最优化计算过程中参数进行迭代更新的条件进行判断；

记录单元，用于记录配准得到的形变场结果、原始图像以及模态变换后的图像。

如图5分别为采用本实施例提供的医学图像配准装置得到的配准后图像以及配准前后叠加图。在本实施例中，当配准过程开始时，需要应用一个全局性的仿射变换以增加配准的准确性。在图像转换单元中初始的仿射变换参数x＝[0，0，0，100，100，0，0]分别代表在x，y方向上的平移、旋转、尺度变换以及裁剪(shearing)。利用配准结果显示模块设置的一个三层的多分辨率方法，同时应用于固定图像和浮动图像，图像归一化单元把图像的灰度值均归一化为[0，1]的区间。在每一个层级的迭代过程中，都应用模态变化的方法。之后，优化方法对已经构建好的能量函数进行优化，能量函数的参数设定为σ₁＝σ₂＝σ＝1。在最优化计算模块中设置最大迭代次数都设定为200次。

在刚体实验中，选取了MRI-T1序列脑部图像和MRI-PD序列脑部图像分别作为固定图像和浮动图像。浮动图像是在固定图像的基础上人为旋转了10°，同时在x和y方向上分别有13mm和17mm的位移。图像大小均为221×257。

实施例3：

本实施例提供一种平台以及一种存储介质。

一种平台包括处理器以及用于存储处理器可执行程序的存储器，所述处理器执行存储器存储的程序时，实现上述实施例1的医学图像配准方法。

一种存储介质，该存储介质为计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时，实现上述实施例1的医学图像配准方法。

本实施例中的存储介质可以是磁盘、光盘、计算机存储器、只读存储器(ROM，Read-Only-Memory)、随机存储器(RAM，Random-Access-Memory)、U盘、移动硬盘等介质。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种基于MRF模型的多模态医学图像配准方法，其特征在于，具体步骤包括：

输入待配准的两幅图像，分别记为固定图和浮动图；

构建马尔可夫随机场模型的能量函数；

将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；

对最终得到的配准结果进行显示。

2.根据权利要求1所述的一种基于MRF模型的多模态医学图像配准方法，其特征在于，所述构建马尔可夫随机场模型的能量函数，包括：

S201、获取输入的待配准图像；

S202、对浮动图进行处理，并将固定图与处理好的浮动图进行归一化；

S203、采用权重图的形式对固定图和浮动图的对应的特殊位置进行标记；

S204、计算固定图和浮动图每个像素之间的距离值，将得到的相似度函数作为数据项；

S205、生成能量函数的平滑项；

S206、根据数据项和平滑项得到马尔可夫随机场模型的能量函数；

S207、将马尔可夫随机场模型的能量函数进行转换；

S208、如果遍历完整图像的所有像素点，则结束；否则转到S204；

S209、根据固定图和处理好的浮动图的归一化结果更新一次联合直方图；

S210、对固定图的一个像素进行迭代，搜索在固定图中重叠次数最多的像素点；

S211、对浮动图的一个像素进行迭代，搜索在浮动图中重叠次数最多的像素点；

S212、判断是否遍历完整图像的所有像素点；

S213、获得模态变换后的固定图与浮动图。

3.根据权利要求2所述的一种基于MRF模型的多模态医学图像配准方法，其特征在于，步骤S202中对浮动图进行处理包括：平移、翻转、旋转以及在x-y方向上的尺度变换；

步骤S204中根据构建的距离函数来表示相似度，表达式为：

E_data＝∑_s∈Vθ_s(x_s)＝1/2σ²∑_s∈V[I(x_s)-J(T(x))]²

其中，I(x)和J(x)分别表示二维空间中的固定图和浮动图，T(·)表示像素点x从固定图I(x)到浮动图J(x)的一个变换模型，在整个定义域Ω中，其中图像I(x)的噪声模型的期望是J(T(x))，σ²表示方差；

给定图像的定义域，将得到的四个图像矩阵结合得到数据项，表达式为：

E_data＝1/2σ₁ ²∑_s∈V[I(x_s)-J(T(x))]²+1/2σ₂ ²∑_s∈V[I(x_s)-J_t(T(x))]²

其中，σ₁和σ₂为高斯白噪声模型的标准差；

在步骤S205中，对于一个领域系统N_i而言，平滑项的相关计算表达式表示为：

E_smooth＝λ∑_s,t∈Niθ_st(x_s,x_t)＝λ∑_s∈V∑_t∈θln[1+(x_s-x_t)²]

其中，λ表示一个权重值，函数θ_st包括分段函数、多项式函数、对数函数；

所述步骤S206中得到的能量函数表示为：

E_mrf(X|θ)＝E_data+E_smooth

在根据步骤S204和S205得到的数据项和平滑项之后，根据马尔可夫随机场模型理论可以得到一个具体的能量函数模型，表示为：

E_mrf＝E_data+E_smooth＝1/2σ₁ ²∑_s∈V[I(x_s)-J(T(x))]²

+1/2σ₂ ²∑_s∈V[I(x_s)-J_t(T(x))]²+λ∑_s,t∈Niln[1+(x_s-x_t)²]

当图像完全配准之后，马尔可夫随机场的能量函数值将会最小；

在步骤S209中，联合直方图更新的公式为：

其中，N表示联合直方图的样条数；x表示当前像素点的位置；表示小于c的最大正整数值；联合直方图H(i,j)是一个二维矩阵；灰度值i在图I(x)以及灰度值j在图像J(x)中相关点具体的数目用H(i,j)来表示；

模态变换公式表示为：

4.根据权利要求1所述的一种基于MRF模型的多模态医学图像配准方法，其特征在于，所述将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数，包括：

S301、初始化sigmoid函数以及优化模块的常数；

S302、初始化步长a₀，根据公式计算出a₁；

S303、根据最终的能量函数，计算固定图和浮动图每个像素及搜索方向；

S304、对固定图和浮动图的像素点分别进行迭代更新；

S305、如果x_i+1＜x_i成立，转到S310；否则转到S306；

S306、如果Δ＝|a_i-a_i-1|≥ε转到S307；否则转到S308；

S307、更新a_i+1＝[a_i·f(-d_i ^T·d_i-1)]⁺；

S308、更新a_i+1＝α·a_i；

S309、令i＝i+1，重新转到S304；

S310、如果遍历完图像所有像素，转到S311；否则转到S304；

S311、记录配准得到的形变场结果、原始图像以及模态变换后的图像，流程结束。

5.根据权利要求1所述的一种基于MRF模型的多模态医学图像配准方法，其特征在于，所述对配准结果进行显示的步骤中，采用一种多分辨率的配准方法显示最终的配准结果，配准过程为：先对一个16×16的子网格像素点进行配准，将得到的形变场和原始图像调整到32×32的子网格再进行配准，以此类推，当配准的网格达到待配准图像的原始大小时配准结束，对配准结果进行显示。

6.一种基于MRF模型的多模态医学图像配准装置，其特征在于，所述装置包括：

图像生成模块，用于输入待配准的两幅图像，分别记为固定图和浮动图；

能量函数构建模块，用于构建马尔可夫随机场模型的能量函数；

最优化计算模块，用于将能量函数输入到改进的优化模型进行最优化计算，得到最终变换参数；

配准结果显示模块，用于对最终得到的配准结果进行显示。

7.根据权利要求6所述的一种基于MRF模型的多模态医学图像配准装置，其特征在于，在所述能量函数构建模块中，包括：

图像转化单元，用于将浮动图进行任意的平移、翻转、旋转以及在x,y方向上的尺度变换；

图像归一化单元，用于将固定图以及浮动图进行归一化处理；

离散标记图生成单元，用于记录固定图和浮动图对应特殊位置的标记，并转化为权重图的形式；

数据项生成单元，用于通过计算固定图和浮动图每个像素之间的距离值，得到相似度函数作为数据项；

平滑项生成单元，用于生产能量函数的平滑项，作为函数模型的约束项；

能量函数生成单元，用于根据平滑项和数据项生产对应的MRF能量函数；

能量函数转化单元，用于将模型转化为最小化能量函数的相关约束；

模态变换单元，用于对固定图和浮动图进行一次基于直方图计算的模态变换；

8.根据权利要求6所述的一种基于MRF模型的多模态医学图像配准装置，其特征在于，在所述最优化计算模块中，包括：

初始化单元，用于初始化sigmoid函数、优化模块的常数以及初始化步长a₀，并根据步长公式计算出a₁；

像素计算单元，用于根据最终的能量函数，计算固定图和浮动图每个像素及搜索方向；

更新单元，用于对固定图和浮动图的像素点以及步长分别进行迭代更新；

判断单元，用于对最优化计算过程中参数进行迭代更新的条件进行判断；

记录单元，用于记录配准得到的形变场结果、原始图像以及模态变换后的图像。

9.一种平台，包括处理器以及用于存储处理器可执行程序的存储器，其特征在于，所述处理器执行存储器存储的程序时，实现权利要求1-5任一项所述的医学图像配准方法。

10.一种存储介质，存储有程序，其特征在于，所述程序被处理器执行时，实现权利要求1-5任一项所述的医学图像配准方法。