CN104050666B - 基于分割的脑部mr图像配准方法 - Google Patents

Abstract

本发明提供一种脑部MR图像配准方法。先将模板脑部MR图像和目标脑部MR图像在经过分割之后，会得到三个明显的组织区域，对于分割出来的组织区域，我们首先计算对应组织之间的相似性，对于不同相似度的组织区域，在进行配准时，我们使用不同的分层等级。如果对应组织区域相似度较高，我们在配准时使用较少的分层以及较高的网格分辨率就可以得到很好的配准结果。当对应组织区域配准完成后，目标图像就可以完整配准到模板图像。

Description

基于分割的脑部MR图像配准方法

技术领域

本发明属于图像处理技术领域，主要基于图像分割的MR脑部图像配准方法。

背景技术

图像配准就是将不同时间、不同成像设备或不同条件下获取的两幅或多幅图像进行匹配、叠加的过程。图像配准技术在遥感领域，医学图像领域，制图领域，和计算机视觉领域都有很高的应用。一般情况下，我们将其主要应用归为四类：多视角分析，多时态分析，多模态分析，场景到模型配准。由于图像多样性以及图像本身不同程度的退化，我们很难找到一种通用的的方法完成所有的配准任务。每一种方法不仅仅需要考虑图像的几何形变，还需要考虑放射形变和噪声干扰等因素。大部分图像配准方法主要有四个步骤：特征检测，特征匹配，变换模型评估，图像重采样和变换。详见文献：BarbaraZitova，JanFlusser.“Imageregistration methods：asurvey”，ImageandVisionComputing21(2003)977-1000所述。

对于脑部核磁共振MR图像配准，需要将目标脑部MR图像与模板脑部MR图像进行匹配、叠加。由于MR图像本身的噪声、灰度不均匀以及弱边界的影响，如何将目标图像完整配准到模板图像是个有待解决的问题。

发明内容

本发明所要解决的技术问题是，提供一种脑部MR图像配准方法。

本发明为解决上述技术问题所采用的技术方案是，基于分割的脑部MR图像配准方法，包括以下步骤：

步骤一，对目标图像I和模板图像I′进行图像分割，分别得到目标图像和模板图像的脑脊液、灰质、白质三个组织区域I_CSF、I_GM、I_WM和I′_CSF、I′_GM、I′_WM；

步骤二，对目标图像的组织区域I_CSF、I_GM、I_WM和模板图像的组织区域I′_CSF、I′_GM、I′_WM分别一一对应采用互信息的方法进行相似测量；

步骤三，对于相似度较高的对应组织区域，设置较低的分层等级参数进行FFD配准，对于相似度较低的对应组织区域，设置较高的分层等级参数进行FFD配准；

步骤四，利用各组织区域对应的分层等级参数，将目标图像组织区域I_CSF、I_GM、I_WM配准到模板图像的组织区域I′_CSF、I′_GM、I′_WM；

步骤五，对于目标图像配准好的组织区域进行组合，对于重叠区域采用融合方法，对于漏空区域采用三线性差值方法得到最终的配准图像。

本发明的有益效果是，先将模板脑部MR图像和目标脑部MR图像在经过分割之后，会得到三个明显的组织区域，对于分割出来的组织区域，我们首先计算对应组织之间的相似性，对于不同相似度的组织区域，在进行配准时，我们使用不同的分层等级。如果对应组织区域相似度较高，我们在配准时使用较少的分层以及较高的网格分辨率就可以得到很好的配准结果。当对应组织区域配准完成后，目标图像就可以完整配准到模板图像。

附图说明

图1为实施例流程图。

具体实施方式

为了方便地描述实施例内容，对现有术语进行介绍：

偏差场。相同组织的不同局部区域的均值，方差有较大的偏差，尤其是当高场强的核磁共振设备出现后，这一问题更加突出。一般采用的偏差场模型为乘性偏差场。

梯度下降法。利用负梯度方向来决定每次迭代的新的搜索方向，使得每次迭代能是优化的目标函数逐步减小。

B-样条。B-样条是样条曲线一种特殊的表示形式。它是B-样条基曲线的线性组合。B-样条是贝兹曲线的一种一般化，可以进一步推广为非均匀有理B样条，使得我们能给更多一般的几何体建造精确的模型。

仿射变换。指在几何中，一个向量空间进行一次线性变换并接上一个平移，变换为另一个向量空间。

互信息。互信息是信息论里一种有用信息的信息度量，它是指两个事件集合之间的相关性，两个事件X和Y的互信息定义为：其中H(X),H(Y)分别代表事件X,Y的熵，且有H(X)＝-∑p(X)logp(X)，H(Y)＝-∑p(Y)logp(Y)H(X,Y)是联合熵，其定义为：H(X,Y)＝-∑p(X,Y)logp(X,Y)，其中p(X,Y)为概率。

本发明是在脑部MR图像分割基础上，应用FFD模型将目标脑部MR图像各个组织区域(白质，灰质，脑脊液)配准到模板脑部MR图像对应的各个组织区域上。然后对目标图像各个组织的配准结果进行综合得到最终的配准结果。具体流程如图1所示：

步骤一，对目标图像和模板图像按照本专利提出的方法先进行图像分割，分别得到目标图像和模板图像的脑脊液、灰质、白质三个组织区域I_CSF、I_GM、I_WM和I'_CSF、I'_GM、I'_WM。

步骤二，对目标图像的组织区域I_CSF、I_GM、I_WM和模板图像的组织区域I'_CSF、I'_GM、I'_WM分别一一对应采用互信息的方法进行相似测量。

步骤三，对于相似度较高的对应组织区域，在进行随后的FFD配准时，我们设置较低的分层等级参数，对于相似度较低的对应组织区域，在进行配准是设置较高的分层等级参数。

步骤四，利用前一等级设置的分层等级参数，将目标图像组织区域I_CSF、I_GM、I_WM配准到模板图像的组织区域I'_CSF、I'_GM、I'_WM。

步骤五，对于目标图像配准好的组织区域进行组合，对于重叠区域采用融合方法，对于漏空区域采用三线性差值方法得到最终的配准图像。

现有的图像分割方法只要是能分割出脑补三组织区域的均可适用于本发明。

常见的图像分割方法有：表面自适应方法，基于直方图的方法，高频最大化方法，信息最小化方法，直方图匹配方法，滤波方法。目前有很多方法用于自动和半自动分割，但是，其中大部分分割效果并不太理想，因为图像本身的噪声，灰度不均匀，以及弱边界等影响。优选的，实施例提供一种先的，特别适用于脑部MR图像分割的方法，如下所示：

本实施例的图像分割方法在现有基于划分的聚类FCM算法(模糊C-均值)的基础上，引入偏差场校正的方法，定义图像中三维像素点x的像素值I(x)如下：

I(x)＝b(x)J(x)+n(x)

其中b(x)为待评估偏差场，其中g₁,…,g_M为正交基函数，即有∫_Ωg_i(x)g_j(x)dx＝δ_ij，其中δ_ij＝0对于i≠j，δ_ij＝1对于i＝j，J(x)为要恢复到的真实图像的三维像素点x的像素值，其中n为加性噪声。

目标函数定义如下：

我们令c＝(c₁,…,c_N)^T，W＝(w₁,…,w_M) ^T，U(x)＝(u₁(x),…,u_N(x))^T，G(x)＝(g₁(x),…g_M(x))^T，然后我们可以将能量函数改写为：

F(U,c,W)＝∫_Ω|I(x)-(W^TG(x)(c^TU(x)))|²dx

T表示转置，b(x)为像素点x的估偏差场，W表示正交基函数系数，G(x)表示正交基函数，M为正交基函数个数，I(x)表示三维像素点x在原始图像中的像素值，J(x)表示三维像素点x在真实图像中的像素值，组织区域类型总数N＝3，Ω表示整个图像区域，c_i表示第i个组织区域的平均灰度值，u_i(x)表示三维像素点x是否属于第i个组织区域的二值标识，0表示不属于，1表示属于；

采用梯度下降法对目标函数做优化处理，相应各个参数的更新如下：

其中 表示取函数最小时i的值；

其中A＝∫_ΩG(x)G(x)^TJ²(x)dx，V＝∫_ΩI(x)G(x)J(x)dx

更新过程如下：

步骤一：初始化c，W和U；

步骤二：更新c到

步骤三：更新W到

步骤三：更新U到

当目标函数更新到给定的最优化条件时(目标函数达到最小值或当前迭代次数达到最大值)，停止迭代。否则转到步骤二继续迭代。

依据当前更新的组织区域的二值标识U完成组织区域的划分，目标图像和模板图像在经过分割之后，会得到三幅图像I_CSF、I_GM、I_WM和I'_CSF、I'_GM、I'_WM，每幅图像对应元目标图像和模板图像的相应组织区域为脑脊液，灰质，白质。

对于图像配准，实施例采用经典的基于B-样条(B-splines)的自由变形算法FFD，该算法使用分层的控制点网格(Φ¹,…,Φ^L)，使其适应不同分辨率等级的需要。算法使用传统的互信息的方法度量目标图像与模板图像的相似性。当算法定义的目标函数达到最优时，就可以得到预期的配准结果。

FFD配准算法描述如下：

图像配准的目的是找到最优的变换T:(x,y,z)→(x',y',z')，我们定义变换T由全局变换T_global和局部变换T_local组成：

T(x,y,z)＝T_global(x,y,z)+T_local(x,y,z)

对全局变换使用现有的仿射变换

对于局部变换，我们使用基于B-样条的FFD模型，FFD通过控制网格点的变化来控制目标的形变。定义Φ为n_x×n_y×n_z的网格，φ_i,j,k为控制节点。FFD的数学公式可以用1-DB-样条的3-D张量积来表示：

其中 表示向下取整，B_l，B_m，B_n代表B-样条的基函数。且有：

B₀(u)＝(1-u)³/6

B₁(u)＝(3u³-6u²+4)/6

B₂(u)＝(-3u³+3u²+3u+1)/6

B₃(u)＝u³/6

采用分层的不同分辨率的网格点(Φ¹,…Φ^L)，每一层的控制网格Φ^l都通过FFD模型得到一个局部变换每一个不同分辨率的局部变换求和得到局部变换

为了使FFD变换在变换过程中更加平滑，我们引入了一个约束项:

其中,表示求导，T为将目标图像匹配到模板图像的变换，λ为预设权值，图像分辨率V＝N_x×N_y×N_z，其中N_x,N_y,N_z分别代表图像在x、y、x三个维度上的大小；该约束项只对局部变换起作用，对全局的放射变换不起作用。

算法使用传统的互信息的方法进行相似性度量：

定义目标函数；

C＝-C_similarity(I′(x),T(I(x)))+λC_smooth

I′(x)表示模板图像中三维像素点x的像素值，T(I(x))表示目标图像中三维像素点x经变换T后的像素值，C_similarity为计算相识度的函数；

配准的具体步骤如下：

步骤一，通过计算相似性度量，找到最优的全局变换参数T_global(x,y,z)。

步骤二，设置分层等级参数L。

步骤三，初始化控制网格点参数Φ^l(n_x×n_y×n_z)，其中l＝1,…,L。

步骤四，确定目标函数，根据梯度下将法求梯度向量Θ代表全局仿射变换参数，在最开始时可直接求出；

步骤五，计算判断是否大于指定阈值ε，如果是则更新||·||表示范数，q为预设调节变量，否则令l＝l+1，如果l＝L，则配准结束，如果l<L，则转步骤三继续迭代。

Claims (3)

1.基于分割的脑部MR图像配准方法，其特征在于，包括以下步骤：

步骤一，对目标图像I和模板图像I'进行图像分割，分别得到目标图像和模板图像的脑脊液、灰质、白质三个组织区域I_CSF、I_GM、I_WM和I'_CSF、I'_GM、I'_WM；

步骤二，对目标图像的组织区域I_CSF、I_GM、I_WM和模板图像中对应的组织区域I'_CSF、I'_GM、I'_WM采用互信息的方法进行相似度测量；

步骤三，根据测量的相似度针对不同组织区域的设置自由变形算法FFD配准的分层等级参数，相似度越高的组织区域设置越低的分层等级参数，根据设置的分层等级参数将目标图像组织区域I_CSF、I_GM、I_WM单独配准到模板图像中对应的组织区域I'_CSF、I'_GM、I'_WM；

步骤四，对于目标图像配准好的组织区域进行组合，对于重叠区域采用融合方法，对于漏空区域采用三线性差值方法得到最终的配准图像。

2.如权利要求1所述基于分割的脑部MR图像配准方法，其特征在于，步骤一的具体方法为：

1)确定目标函数F(U,c,W)＝∫_Ω|I(x)-(W^TG(x)(c^TU(x)))|²dx，

其中，I(x)表示三维像素点x在原始图像中的像素值，I(x)＝b(x)J(x)+n(x)，b(x)为像素点x的待评估偏差场，J(x)表示三维像素点x在真实图像中的像素值，n(x)为加性噪声；c＝(c₁,…,c_N)^T，U(x)＝(u₁(x),…,u_N(x))^T，W＝(w₁,…,w_M)^T，G(x)＝(g₁(x),…g_M(x))^T，T表示转置，W表示正交基函数系数，G(x)表示正交基函数，M为正交基函数个数，N为组织区域类型总数N＝3，Ω表示整个图像区域；

2)采用梯度下降法对目标函数求极小值，各参数的更新如下：

2-1)初始化图像中组织区域的二值标识U、正交基函数系数W、组织区域的平均灰度值c；

2-2)更新图像中组织区域的二值标识U中各元素值，i＝1…,N，

其中， 表示取函数最小时i的值；

2-3)更新各组织区域的平均灰度值：

i＝1…,N，c_i表示第i个组织区域的平均灰度值，u_i(x)表示三维像素点x是否属于第i个组织区域的二值标识，0表示不属于，1表示属于；

2-4)更新正交基函数系数，W＝A^-1V，其中，中间量A＝∫_ΩG(x)G(x)^TJ²(x)dx，中间量V＝∫_ΩI(x)G(x)J(x)dx；

2-5)将当前更新的组织区域的二值标识U、正交基函数系数W、组织区域的平均灰度值c代入目标函数，判断是否目标函数达到最小值或当前迭代次数达到最大值，如是则停止迭代，依据当前更新的组织区域的二值标识U完成组织区域的划分，否则转到步骤2-2)继续迭代更新各参数。

3.如权利要求1所述基于分割的脑部MR图像配准方法，其特征在于，步骤三的具体方法为：

1)根据测量的相似度设置分层等级参数L；

2)用分层等级参数L初始化不同分辨率的网格点Φ^l＝Φ¹,…Φ^L，其中l∈(1,…,L)，初始化l＝1；

3)确定目标函数C＝-C_similarity(I′(x),T(I(x)))+λC_smooth作为目标函数，根据梯度下降法求梯度向量其中，C_smooth(T)为约束项， 表示求导，T为将目标图像匹配到模板图像的变换，λ为预设权值，图像分辨率V＝N_x×N_y×N_z，其中N_x,N_y,N_z分别代表图像在x、y、z三个维度上的大小；I′(x)表示模板图像中三维像素点x的像素值，T(I(x))表示目标图像中三维像素点x经变换T后的像素值，C_similarity为计算相似度的函数，Θ表示全局仿射变换参数；

4)判断梯度向量是否大于指定阈值ε，如是，则更新||·||表示范数，q为预设调节变量，否则，令l＝l+1，如果l＝L，则配准结束，将当前变换T作为目标图像匹配到模板图像的最优变换，如果l<L，则转步骤3)继续迭代。