一种基于自适应邻域选择的MR图像配准方法与流程

本发明涉及图像处理技术领域，特别是涉及一种mr图像配准方法。

背景技术

正常情况下，配准后的解剖结构图像与参考图像应具有相同的拓扑结构，不会产生新的组织，原有的结构不会消失，形变域也不会出现撕裂，折叠，空洞等不合理的物理结构。微分同胚是一种平滑，连续和可逆的变换，能够使解剖结构图像变形后的拓扑结构保持不变，在医学图像配准中是一个重要的应用，如下。

在文献[1](marslands,twiningcj.constructingdiffeomorphicrepresentationsforthegroupwiseanalysisofnonrigidregistrationsofmedicalimages[j].ieeetransactionsonmedicalimaging,2004,23(8):1006-1020.)，marsland等把微分同胚变换看作时变速度场，用测地插值样条基构造变换。

在文献[2](begmf,millermi,trouvéa,etal.computinglargedeformationmetricmappingsviageodesicflowsofdiffeomorphisms[j].internationaljournalofcomputervision,2005,61(2):139-157.)，beg等提出lddmm算法在速度场空间利用基于欧拉-拉格朗日方程的变分方法求解大形变的微分同胚变换，采用梯度下降法更新形变参数。

在文献[3](ashburnerj,fristonkj.diffeomorphicregistrationusinggeodesicshootingandgauss–newtonoptimisation[j].neuroimage,2011,55(3):954-967.)，ashburner等将[2]做了改进，任意时刻的速度场都用初始速度场来表示，配准的每次迭代都利用初始形变值来计算，减少了对内存和外存的要求。

在文献[4](ashburnerj.afastdiffeomorphicimageregistrationalgorithm[j].neuroimage,2007,38(1):95-113.)，ashburner等提出dartel非时变速度场的配准方法，同胚变换构成了复合运算下的李群，将大形变同胚变换分解为一系列小形变来处理，使可逆变换相对容易，并降低计算代价。

在文献[5](janssensg,jacquesl,dexivryjo,etal.diffeomorphicregistrationofimageswithvariablecontrastenhancement[j].journalofbiomedicalimaging,2011,2011:3.)，janssens等在微分同胚方法中利用累加位移场的可逆性配准对比增强度不同的图像。

在文献[6](arsignyv,commowicko,pennecx,etal.alog-euclideanframeworkforstatisticsondiffeomorphisms[c]//internationalconferenceonmedicalimagecomputingandcomputer-assistedintervention.springerberlinheidelberg,2006:924-931.)，arsigy等借助李群理论提出log-euclidean框架，在微分同胚空间里可以对向量场进行分析统计并且能保持变换的可逆性。

在文献[7](vercauterent,pennecx,perchanta,etal.diffeomorphicdemons:efficientnon-parametricimageregistration[j].neuroimage,2009,45(1):s61-s72.)，vercauteren等基于demons算法，提出一种非参数diffeomorphicdemons配准算法(dd-np)，利用李群理论在连续域上进行空间变换，使得空间变换具有微分同胚性，从而使形变场具有拓扑保持性。

在文献[8](dt-refind:diffusiontensorregistrationwithexactfinite-straindifferential[j].medicalimagingieeetransactionson,2009,28(12):1914-1928.)，yeo等利用有限应变微分对张量方向重定向，把[7]从标量图像配准扩展到了张量图像配准。

如上述，现有基于微分同胚的配准技术都没有考虑在高维微分同胚变换空间中数据的非线性结构，而高维数据通常包含丰富的结构信息，这种结构信息对空间变换的拓扑保持性具有显著影响，最终影响配准结果的准确性。

技术实现要素：

本发明的目的在于提供一种基于自适应邻域选择的mr图像配准方法，对原始的微分同胚demons技术进行改进，引入流形学习方法，提出了一种自适应邻域选择的mr图像配准技术，使得变换空间具有拓扑保持性。

如果局部切空间的线性化程度越高,则高维数据流形的非线性结构信息就越丰富，在空间变换中就可以最大限度上保留流形的局部非线性结构。那么如何能更好地逼近切空间，使得切空间的线性化程度提高呢？由于局部切空间的逼近与样本点的邻域有密切关系，通过流形学习方法自适应地选择邻域大小，形成线性子空间，达到更好地逼近切空间的目的，进而在空间变换中保持图像的拓扑结构，改善这个技术难点引起的问题，最终提高图像配准的精度。

为了实现上述目的，本发明提供了以下技术方案：

本发明提供的一种自适应邻域选择技术，包括以下步骤，

1.1、基于李群和李群的单位元id处的切空间之间的映射关系，获取在局部坐标系中图像任意样本点的形变场u，

设李群的单位元id处的切空间是tg，其d维正交基矩阵为t＝[τ1,τ2…τd]，τi∈r^m,i＝1,2,…,d，在局部坐标系中，任意样本点的形变场u可以表示为下述公式(f-1-1)的线性组合，

其中，uⁱ是给定坐标下形变场u的分量，是线性组合系数；

1.2、采用线性拟合来近似计算邻域点εj的线性逼近；

1.3、通过kernelpca方法，找到一个隐式映射函数φ(si)，把si投影到具有更好性质的特征空间,再在特征空间中寻找可以保持最大差异性的d个低维投影方向[τ1,τ2…τd]，最终得到邻域点εj在局部邻域上的投影u^j，也就是邻域点的坐标；

1.4、基于邻域点εj的坐标u^j计算对应正交基矩阵的局部坐标系u，并用svd分解得到左右奇异值；

1.5、利用左右奇异值计算比值，得到邻域选择标准r。

本发明提供的一种基于自适应邻域选择的mr图像配准方法，包括以下步骤，

2.01、采集两幅mri图像，一幅为参考图像ir，另一幅为浮动图像if；

2.02、从两幅图像的像素中分别提取对称正定的图像特征，构造出高维李群流形结构；

2.03、预先设定一个最小邻域kmin，用k-最近邻法确定一个初始邻域大小为k的单位元邻域集eε^k；

2.04、进行自适应邻域选择，先用svd分解得到奇异值，再利用奇异值计算比值得到邻域选择标准r；

2.05、设定邻域选择标准r的阈值η∈[0,1]，将r与η进行比对，选择邻域大小；若是r≥η且k＞kmin，则移除距离邻域均值最远的点，使得k←k-1，则邻域尺寸缩小为eε^k←eε^k-1，然后返回程序2.4；若是r＜η，则邻域大小为k，然后执行下一程序；

2.06、计算获取邻域集eε^k的扩张邻域尺寸eε^k←eε^k+1；

2.07、计算李群单位元处邻域所对应的切空间的正交基矩阵t，获得形变场u；

2.08、利用形变场u计算参考图像ir中的任意样本点si的雅可比行列式

2.09、通过最小化代价函数，计算得到最优空间变换其中，最小化代价函数的约束条件是雅可比行列式的值且的值恒为正，使得参考图像ir和浮动图像if中的点一一对应，保持拓扑结构；

2.10、通过最优空间变换对浮动图像if进行配准，输出配准图像。

与现有技术相比，本发明具有以下优点：

本发明在李群流形下，用自适应的邻域选择方法去逼近单位元处的切空间，这样能更精确地逼近流形的局部几何结构，以保证空间变换的一一映射性，配准前后的浮动图像仍保持邻接关系，使配准后的图像结构保持拓扑性，最终达到改善精度的目的。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明的mri图像配准框架示意图；

图2为本发明的李群和切空间之间的映射关系图，其中，gl(m)表示m×m矩阵构成李群，gl(m)表示m×m矩阵构成的李代数(线性向量空间)，log()表示对数映射，exp()表示指数映射，si∈gl(m),si∈gl(m)是指数映射使m维流形局部同胚于m维向量空间；

图3为本发明的临床数据的配准结果对比图，其中，w1是白质参考图像，w2是白质浮动图像，w3是dd-np处理白质图像，w4是本发明方法处理白质图像，g1是灰质参考图像，g2是灰质浮动图像，g3是dd-np处理灰质图像，g4是本发明方法处理灰质图像；

图4为本发明的临床数据中大脑白质上rmse的对比图；

图5为本发明的临床数据中大脑白质上rmse的对比图。

具体实施方式

以下结合附图对本发明的实施例进行说明，应当理解，此处所描述的实施例仅用于说明和解释本发明，并不用于限定本发明。

实施例

流形上任一点附近的局部线性结构都能被流形在该点处的切空间描述，切空间是该局部线性化子空间，线性化程度越高，包含的流形的非线性结构越多。

参见图1，为了得到高维空间的结构信息，首先，构造正定对称的协方差矩阵，协方差的空间结构是一个高维黎曼流形，在乘法和可逆条件下形成李群；其次，基于李群可由单位元的一个任意小领域生成，我们利用样本点邻域的局部切空间来表示李群流形的非线性几何结构；再次，通过流形学习方法自适应地选择邻域大小，提高局部切空间的线性化程度越高；最后，最小化代价函数得到最优变换空间，进而配准输出。

综上，本发明要到达的第一个目的是对邻域进行自适应选择；在第一个目的的基础上，第二个目的是进行高精度的mr图像配准。具体方案如下：

本发明提供的一种自适应邻域选择技术，包括以下步骤，

步骤1.1、基于李群和李群的单位元id处的切空间之间的映射关系，参见图2，获取在局部坐标系中图像任意(随机)样本点的切向量空间(形变场)u。

设李群的单位元id处的切空间是tg，其d维正交基矩阵为t＝[τ1,τ2…τd]，τi∈r^m,i＝1,2,…,d，在局部坐标系中，任意样本点的形变场u可以表示为下述公式(f-1-1)的线性组合，

其中，uⁱ是给定坐标下形变场u的分量，是线性组合系数。

步骤1.2、采用线性拟合来近似计算邻域点εj的线性逼近。具体内容如下：

1.2.1、设eε是单位元id的邻域，且eε由k个最近邻点[ε1,ε2,…εk]组成，εj∈sym⁺(m),j＝1,2,…k，sym⁺(m)为一个由m维的正定对称矩阵构成的空间；

1.2.2、邻域点εj的线性逼近用线性拟合来近似，通过下述公式(f-1-2)计算，

其中，t∈r^m×d是构成切空间的正交基矩阵，是k个最近邻的均值。

步骤1.3、通过kernelpca方法，找到一个隐式映射函数φ(si)，把si投影到具有更好性质的特征空间,再在特征空间中寻找可以保持最大差异性(方差)的d个低维投影方向[τ1,τ2…τd]，最终得到邻域点εj在局部邻域上的投影u^j，也就是邻域点的坐标。邻域点εj的坐标u^j通过下述公式(f-1-3)计算：

步骤1.4、基于邻域点εj的坐标u^j计算对应正交基矩阵的局部坐标系u，并用svd分解得到左右奇异值。具体内容如下：

1.4.1、设u＝[u¹,u²,…u^k],u∈r^d×k，是对应正交基矩阵的局部坐标系，根据u^j将其转成带有的计算式，见下述公式(f-1-4)，

1.4.2、用svd(奇异值分解)方法分解得到k个奇异值，按序排列为σ1≥…≥σd≥…σk，于是公式(f-1-4)可以表达为下述公式(f-1-5)，

u＝diag(σ1,σ2,…σd)v^τ(f-1-5)

其中，v^τ是的d个最大奇异值对应的右奇异向量，则t是左奇异向量构成的正交基矩阵。

步骤1.5、利用左右奇异值计算比值，得到邻域选择标准r。具体内容如下：

1.5.1、根据左右奇异值得到邻域选择标准r的比值参数的计算式，见下述公式(f-1-6)和公式(f-1-7)，

以及，

1.5.2、根据比值参数计算比值，即公式(f-1-6)和公式(f-1-7)，构造出邻域选择标准r，见下述公式(f-1-8)，

如上述，通过邻域选择标准实现自适应邻域选取。

下述将对mr图像配准方法进行说明，其是基于如上述的自适应邻域选择技术的mr图像配准技术，包括以下步骤，

步骤2.01、通过核磁共振仪，采集两幅mri图像，一幅为参考图像ir，另一幅为浮动图像if。

步骤2.02、从两幅图像的像素中分别提取对称正定的图像特征，构造出高维李群流形结构。

步骤2.03、预先设定一个最小邻域kmin，用k-最近邻法确定一个初始邻域大小为k的单位元邻域集eε^k。

步骤2.04、进行自适应邻域选择，先用svd分解得到奇异值，再利用奇异值计算比值得到邻域选择标准r。过程如下：

2.04.1、基于李群和李群的单位元id处的切空间之间的映射关系，获取在局部坐标系中任意样本点的形变场u；

2.04.2、采用线性拟合来近似计算邻域点εj的线性逼近；

2.04.3、通过kernelpca方法，找到一个隐式映射函数φ(si)，把si投影到具有更好性质的特征空间,再在特征空间中寻找可以保持最大差异性(方差)的d个低维投影方向[τ1,τ2…τd]，最终得到邻域点εj在局部邻域上的投影u^j，也就是邻域点的坐标；

2.04.4、基于邻域点εj的坐标u^j计算对应正交基矩阵的局部坐标系u，并用svd分解得到左右奇异值；

2.04.5、利用左右奇异值计算比值，得到邻域选择标准r。具体内容参见上述自适应邻域选择技术。

步骤2.05、设定邻域选择标准r的阈值η∈[0,1]，将r与η进行比对，选择邻域大小；若是r≥η且k＞kmin，则移除距离邻域均值最远的点，使得k←k-1，则邻域尺寸缩小为eε^k←eε^k-1，然后返回程序2.4；若是r＜η，则邻域大小为k，然后执行下一程序。

步骤2.06、计算获取邻域集eε^k的扩张邻域尺寸eε^k←eε^k+1。具体内容如下：

计算其中，εj，k+1＜j≤kmax，即初始邻域中没有选中的点；如果满足约束则将εj加回eε，使k←k-1，eε^k←eε^k+1。

步骤2.07、计算李群单位元处邻域所对应的切空间的正交基矩阵t，获得形变场u。

步骤2.08、利用形变场u计算参考图像ir中的任意(随机)样本点si的雅可比行列式

si的雅可比行列式见下述公式(f-2-1)：

步骤2.09、通过最小化代价函数，计算得到最优空间变换其中，最小化代价函数的约束条件是雅可比行列式的值且的值恒为正，使得参考图像ir和浮动图像if中的点一一对应，保持拓扑结构，最终使得精度提高。

最小化代价函数见下述公式(f-2-2)：

步骤2.10、通过最优空间变换对浮动图像if进行配准，输出配准图像。

如上述，找到最优变换空间后，对浮动图像if进行配准，输出配准图像，实现mr图像配准，参见图1-2。在李群流形下，用自适应的邻域选择方法去逼近单位元处的切空间；邻域大小与样本点的曲率密切相关，流形上曲率是高度变化的，曲率大的样本点的邻域应该相对比较小,而流形上曲率小的样本点处的邻域比较大；这样能更精确地逼近流形的局部几何结构，以保证空间变换的一一映射性，配准前后的浮动图像仍保持邻接关系，使配准后的图像结构保持拓扑性，最终达到改善精度的目的。

下述将结合具体实例的图像处理过程进行说明。

步骤1参考图像和浮动图像的数据采集来自西门子核磁共振仪，回波时间te＝2.98ms，回波链长度etl＝1，反转时间ti＝1100ms，重复时间tr＝2530ms，视野fov＝87.5，图像大小448×512×192，20个健康主体大脑图像。

步骤2随机抽选4个主体的冠状面和横断面图像，然后设为两组：冠状面组用于分割大脑白质，横断面组用于分割大脑灰质。每组选一张作为参考图像，另一张为浮动图像。

步骤3构造李群流形以获取高维结构信息：对选取的每个样本像素构造一个136×136的正定对称的协方差矩阵作为图像特征。

步骤3.1提取145维局部特征fi，包括128维的sift描述子，坐标位置(x，y)，灰度值ii(xy),灰度值的一阶梯度，一阶梯度的模和二阶梯度，此外，把样本图像分割为若干个像素的单元(cell)，把梯度方向平均划分为9个区间(bin)，在每个单元里面对所有像素的梯度方向在各个方向区间进行直方图统计，可得到一个9维的特征向量。

步骤3.2通过fi和fi的转置做外积计算，得到如下的正定对称的协方差矩阵si＝fifi^τ。

步骤4对得到的协方差矩阵所形成的高维特征做pca降维处理，减少噪声影响和计算代价。

步骤5降维后的数据根据如前述的邻域自适应选择技术进行自适应邻域选取，获得邻域选取标准r。

步骤6根据如前述的图像配准技术，获得最优变换空间

步骤7对浮动图像进行配准，输出配准图像。

步骤8检验。两组实验各进行20次，最后计算每20次实验的平均均方根误差rmse。

下述将结合临床数据进行说明。

参见图3，临床数据的配准结果对比，第一行(w1-w4)和第二行(g1-g4)分别显示了临床数据的大脑白质和灰质配准在两种算法上的对比可视化结果，可以看出本发明方法比dd-np方法在白质和灰质成像上均更接近参考图像。

参见图4，临床数据中大脑白质上rmse的对比，图4显示在m＝135，k＝14时的本发明方法与dd-np算法在白质图像上的误差与迭代关系。参见图5，临床数据中大脑白质上rmse的对比，图5显示在m＝125，k＝18时的本发明方法与dd-np方法在灰质图像上的误差与迭代关系。可以看出，在第80次迭代以后，本发明方法的误差均小于dd-np方法。

应当理解，本发明上述实施例及实例，是出于说明和解释目的，并非因此限制本发明的范围。本发明的范围由权利要求项定义，而不是由上述实施例及实例定义。