CN103761750A - 一种心肌质点运动图像与心肌纤维走向图像配准方法 - Google Patents

Abstract

本发明提供一种心肌质点运动图像与心肌纤维走向图像配准方法：根据Tagging MRI图像序列获得形变位移场图像和邻域伸缩运动图像及其相对运动方向中心点，并根据原始心脏DTI图像获取心肌纤维走向图像及其纤维方向中心点，对齐两个中心点并确定相对平移量；将两个图像均分多个扇形，获取两图像间相似度最大值对应的相对旋转角度；再在该相对旋转角度下将其中一幅图像相对另一幅图像按不同倍数缩放，求取两图像间最大相似度的放大倍数，根据上述相对平移量、相对旋转角度和相对放大倍数对上述图像变换处理，得到配准结果。本发明将纤维的方向和心肌自主运动的方向作为一对待配准异元方向特征，提供一种在纤维层次上配准的方法。

Description

一种心肌质点运动图像与心肌纤维走向图像配准方法

技术领域

本发明属于医学图像处理技术领域，更具体地，涉及一种心肌质点运动图像与心肌纤维走向图像配准方法。

背景技术

伴随着我国经济的快速增长，人民生活方式的改变，心脏病已成为城乡人口的主要死因之一，它严重威胁着人民的生命健康，给社会带来了巨大的经济压力。深入研究心脏生理病理的本质规律，特别是实现心脏病的早期无创诊断具有十分重大的意义。基于计算机科学及生物医学等领域的新技术对心脏进行建模仿真是目前进行心脏研究的有效手段，它能够借助计算机强大的计算和可视化能力解决复杂的心脏运动模拟问题，从而辅助医生进行心脏的早期诊断，大大降低心脏病患者的死亡率和发病率。

在对心脏进行研究时，往往需要参考多种模态的图像(来源于不同成像模式的图像)所提供的心脏结构、运动及功能等信息，这种情况下一般要进行多模态图像的配准处理，即对其中一幅图像进行一定的几何变换映射到另一幅图像中，使得两幅图像中的相关点达到空间上的一致。对多种模态的心脏影像进行配准，融合诊断心脏疾病所需的图像信息，具有重要的理论意义和广阔的应用前景。

常用来分析心脏功能的医学成像技术包括Tagging MRI(TaggingMagnetic Resonance Imaging,带标记的核磁共振图像)、DTI(DiffusionTensor Imaging，弥散张量成像)、CT(Computerize Tomography，计算机断层扫描)、US(Ultrasound，超声波)、PET(Positron Emission Tomography，正电子发射断层扫描)等等，每种模态的影像反映了复杂心脏运动过程不同侧面的信息。对多种模态的心脏图像进行配准，特别是将能够获取心肌纤维结构的离体DTI图像，与能够提供活体心脏运动信息的Tagging MRI图像进行配准，是建立心肌纤维的运动模型及对心脏功能进行分析研究的有效途径之一。

发明内容

本发明的目的在于提供一种心脏心肌质点运动图像与心脏心肌纤维走向图像在纤维层次上的配准方法，该方法能够得到心脏DTI图像与心脏Tagging MRI在纤维层次上的配准结果，为建立心脏纤维运动力学模型、在纤维层次上评估心脏运动功能打下基础，进而为分析心脏功能、心脏病病理学研究提供客观可靠的依据。

为实现上述目的，本发明提供了一种心肌质点运动图像与心肌纤维走向图像配准方法，该方法包括下述步骤：

第1步，根据Tagging MRI图像序列获得心脏心肌质点的形变位移场图像M_u，也即为心肌的质点运动图像；

第2步，根据上述形变位移场图像M_u计算邻域伸缩运动图像M_f；

第3步，计算上述邻域伸缩运动图像M_f的运动中心点S₁；

第4步，获得原始心脏DTI图像的心肌纤维走向图像M_a；

第5步，计算上述纤维走向图像M_a的中心点S₂；

第6步，将上述邻域伸缩运动图像M_f的运动中心点S₁与上述纤维走向图像M_a的中心点S₂对齐，并确定平移量；

第7步，以对齐后的中心点为圆心向外引多条射线，将上述邻域伸缩运动图像M_f和纤维走向图像M_a分成大小相同的多个扇形；

第8步，计算邻域伸缩运动图像M_f和纤维走向图像M_a之间的相似度度量函数值f，

其中符号E{}表示求期望值，

a_i表示邻域伸缩运动图像M_f中第i个扇形区域内所有像素点形变方向投影到该扇形平分线的垂线上的投影值之和,b_i表示纤维走向图像M_a中第i个扇形区域内所有像素点形变方向投影到该扇形平分线的垂线上的投影值之和,n为邻域伸缩运动图像M_f和纤维走向图像M_a中扇形的个数；

第9步，将邻域伸缩运动图像M_f和纤维走向图像M_a相对旋转一个扇形角度大小，并按第8步的方法计算旋转后邻域伸缩运动图像M_f和纤维走向图像M_a之间的相似度度量函数值f，直到邻域伸缩运动图像M_f和纤维走向图像M_a之间相对旋转一周；

第10步，当邻域伸缩运动图像M_f和纤维走向图像M_a之间相对旋转一周得到每个旋转角度下的相似度度量函数值f时，确定上述相似度度量函数值f中的最大值，以及上述相似度度量函数值f的最大值所对应的旋转角度α；

第11步，相对邻域伸缩运动图像M_f和纤维走向图像M_a之间的初始位置，将邻域伸缩运动图像M_f和纤维走向图像M_a相对旋转上述角度α，并将旋转后的邻域伸缩运动图像M_f和纤维走向图像M_a中的任一图像相对另一图像以中心点为中心依次放大K₁,K₂,…,K_m倍，计算每次放大后邻域伸缩运动图像M_f和纤维走向图像M_a之间的相似度度量函数值g，得到m个相似度度量函数值，并确定这m个相似度度量函数值中的最大值g_t，以及最大值g_t所对应的放大系数；其中，m为放大的次数，相似度度量函数值g具体计算表达式为：

$g=\underset{i=1}{\overset{t}{\mathrm{\Σ}}}\left|\stackrel{\→}{u_i}\right|.\left|\stackrel{\→}{v_i}\right|.\cos (\stackrel{\→}{u_i},\stackrel{\→}{v_i})$

其中t为邻域伸缩运动图像M_f和纤维走向图像M_a重合的点的个数，与

分别为两幅图像中相同位置的两像素点对应的方向向量，

与

分别为对应方向向量的模，

表示求两向量与

的夹角的余弦值；

第12步，配准搜索计算本次迭代结束，判断本次最大相似度度量函数值g_t与上次最大相似度度量函数g_t-1之差的绝对值是否小于相似度阈值r₁，同时判断扇形的角度是否小于角度阈值r₂，若上述两个条件之一成立则配准搜索计算过程结束，同时转向第14步，否则转向第13步继续搜索；

第13步，减小邻域伸缩运动图像M_f和纤维走向图像M_a对应的等分扇形角度的大小，并转向第7步；

第14步，根据前述步骤中确定的图像间平移量，以及最佳旋转角度和最优放大系数，将邻域伸缩运动图像M_f或纤维走向图像M_a按照这些参数进行转换，得到配准结果。

进一步地，所述第1步具体为：根据Tagging MRI图像序列计算获得心脏心肌质点形变位移场图像也即心脏心肌质点运动图像，心肌在拉伸与收缩的过程中，像素点(p,q)的形变用二维向量

来表示，其中x，y分别为该像素点在水平与垂直方向上的位移量，该向量被定义为形变向量，心脏Tagging MRI图像中所有像素的形变向量构成了该图像的形变位移场，记整个心脏Tagging MRI图像的形变位移场图像为M_u。

进一步地，所述第2步具体包括：

Step1：若Tagging MRI图像所对应的时刻处于心脏收缩期，则转向step2；否则，若Tagging MRI图像所对应的时刻处于心脏舒张期，则转向step3；

Step2：根据每个像素(p,q)的邻域相对形变场R_pq，按照公式

$\stackrel{\→}{L_c}=\underset{\stackrel{\→}{L}}{\arg \max }\left(\frac{1}{2}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\right|\stackrel{\→}{r_k}|.\cos {\mathrm{\θ}}_{(\stackrel{\→}{r_k},\stackrel{\→}{L})}.[1+\frac{\left|\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|\right|}{\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|}\left]\right)$

计算该像素的邻域收缩主轴方向向量

获得图像的邻域收缩主轴方向场M_c；转Step4；其中定义为：心脏Tagging MRI图像中像素点(p,q)的邻域相对形变场中，将相对形变向量

(k=1,2,3,…K)向过中心像素任意方向的直线L投影，方向向量记为

其中朝向中心像素的投影向量长度之和最大的直线方向，被称为该像素的邻域收缩主轴方向，其方向向量记为

K为邻域像素点个数；邻域收缩主轴方向场M_c为：所有像素的邻域收缩主轴方向向量组成了该图像的邻域收缩主轴方向场；

为向量

与向量

的夹角；

和

分别为中心像素点到点A^(k)、B^(k)的矢量，点A^(k)、B^(k)分别为第k个邻域像素的相对形变向量

投影到过原点的直线L上的投影向量的起点与终点；||符号表示求模；

Step3：根据每个像素(p,q)的邻域相对形变场R_pq，按照公式

$\stackrel{\→}{L_s}=\underset{\stackrel{\→}{L}}{\arg \max }\left(\frac{1}{2}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}\right|\stackrel{\→}{r_k}|.\cos {\mathrm{\θ}}_{(\stackrel{\→}{r_k},\stackrel{\→}{L})}.[1-\frac{\left|\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|\right|}{\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|}\left]\right)$

计算该像素的邻域拉伸主轴方向向量

获得图像的邻域拉伸主轴方向场M_s；转Step4；其中

定义为：心脏Tagging MRI图像中像素点(p,q)的邻域相对形变场中，将相对形变向量

(k=1,2,3,…K)向过中心像素任意方向的直线L投影，方向向量记为其中远离中心像素的投影向量长度之和最大的直线方向，被称为该像素的邻域拉伸主轴方向，其方向向量记为

K为邻域像素点个数；邻域拉伸主轴方向场M_s为：所有像素的邻域拉伸主轴方向向量组成了该图像的邻域拉伸主轴方向场；

为向量

与向量

的夹角；

和

分别为中心像素点到点A^(k)、B^(k)的矢量，点A^(k)、B^(k)分别为第k个邻域像素的相对形变向量

投影到过原点的直线L上的投影向量的起点与终点；||符号表示求模；

Step4：若心肌在该时刻处于收缩期内，则该Tagging MRI图像所对应的邻域伸缩运动图像M_f＝M_c；反之，若心肌在该时刻处于舒张期内，则M_f＝M_s。

其中邻域相对形变场R_pq的计算方式如下：将心脏Tagging MRI图像中像素点(p,q)的第k个邻域像素记为(p,q)^(k)，该邻域点相对于像素点(p,q)的相对形变向量为两像素点形变向量的差，记为所有邻域像素的相对形变向量组成该像素的邻域相对形变场，记为R_pq。

进一步地，所述第3步具体包括：过图像M_f上的每一像素点作一条垂直于该像素点形变方向的直线，计算运动中心点S₁到该直线的距离，同时将该距离乘以相应像素点形变大小的幅值，将所有像素点计算得到的该值相加，使得该值最大化的点即为所求的运动中心点，记图像M_f的运动中心点为S₁。

本发明提供的配准方法是一种心脏多模态图像配准方法，能够解决心脏DTI图像与心脏Tagging MRI图像在纤维层次上的配准问题，DTI图像能够反映离体心脏的纤维结构，但不能反映活体心脏的运动信息，通过将心脏DTI图像与能提供活体心脏运动信息的Tagging MRI图像进行配准，来估计心肌纤维的运动信息。本发明方法根据心肌纤维具有沿着纤维伸展方向进行自主收缩与舒张运动的特点，将纤维的方向和心肌自主运动的方向作为物理含义截然不同而又密切相关的一对待配准的特征，能够提供一种将DTI图像和Tagging MRI图像在纤维层次上配准的方法，为建立心脏纤维运动力学模型、在纤维层次上评估心脏运动功能打下基础，进而为分析心脏功能、心脏病病理学研究提供客观可靠的依据。

附图说明

图1为本发明心肌质点运动图像与心肌纤维走向图像配准方法的流程图；

图2为本发明形变位移场示意图；

图3为本发明邻域相对形变场示意图；

图4为本发明邻域伸缩运动图像示意图；

图5为本发明计算运动中心点的示意图；

图6为本发明纤维走向图像示意图；

图7为本发明图像旋转示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明心肌质点运动图像与心肌纤维走向图像配准方法，包括以下步骤：

（1）根据Tagging MRI图像序列计算获得心脏心肌质点形变位移场图像也即心脏心肌质点运动图像。心肌在拉伸与收缩的过程中，像素点(p,q)的形变用二维向量

来表示，其中x，y分别为该像素点在水平与垂直方向上的位移量。该向量被定义为形变向量，心脏Tagging MRI图像中所有像素的形变向量构成了该图像的形变位移场，记该形变位移场为M_u，获取的形变位移场图像如图2所示。

（2）获取心肌质点形变位移场图像的邻域伸缩运动图像。将有关定义描述如下：

定义1心脏Tagging MRI图像中像素点(p,q)的第k个邻域像素记为(p,q)^(k)，该邻域点相对于像素点(p,q)的相对形变向量为两像素点形变向量的差，记为

所有邻域像素的相对形变向量组成该像素的邻域相对形变场，记为R_pq。

分析某像素的邻域相对形变场，以中心像素点为平面坐标原点O建立平面坐标系，如图3所示(以8邻域为例)。将第k个邻域像素的相对形变向量

投影到过原点的直线L上，投影向量的起点与终点分别记为A^(k)、B^(k)。根据投影始点、终点相对于中心像素点的位置，可判断该点相对形变是朝向还是远离中心像素点。即投影始点A^(k)到原点的距离若大于投影终点B^(k)到原点的距离，则

的投影向量朝向中心像素点；反之，投影始点A^(k)到原点的距离若小于投影终点B^(k)到原点的距离，则

的投影向量远离中心像素。本实施例中以8邻域为例，当然邻域的个数也可以为其他的值。

定义2心脏Tagging MRI图像中像素点(p,q)的邻域相对形变场中，将相对形变向量(k=1,2,3,…8)向过中心像素任意方向的直线L投影(方向向量记为

)，其中朝向中心像素的投影向量长度之和最大的直线方向，被称为该像素的邻域收缩主轴方向，其方向向量记为

可以由公式(1)计算。所有像素的邻域收缩主轴方向向量组成了该图像的邻域收缩主轴方向场，记为M_c

$\stackrel{\→}{L_c}=\underset{\stackrel{\→}{L}}{\arg \max }\left(\frac{1}{2}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}\right|\stackrel{\→}{r_k}|.\cos {\mathrm{\θ}}_{(\stackrel{\→}{r_k},\stackrel{\→}{L})}.[1+\frac{\left|\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|\right|}{\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|}\left]\right)---\left(1\right)$

其中，

为向量

与向量

的夹角，||符号表示求模。

定义3像素点(p,q)的邻域相对形变场中，将相对形变向量

(k=1,2,3,…8)向过中心像素任意方向的直线L投影(方向记为

)，其中远离中心像素的投影向量长度之和最大的直线方向，被称为该像素的邻域拉伸主轴方向，其方向向量记为

可由公式2计算。所有像素的邻域拉伸主轴方向向量组成了该图像的邻域拉伸主轴方向场，记为M_s。

$\stackrel{\→}{L_s}=\underset{\stackrel{\→}{L}}{\arg \max }\left(\frac{1}{2}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}\right|\stackrel{\→}{r_k}|.\cos {\mathrm{\θ}}_{(\stackrel{\→}{r_k},\stackrel{\→}{L})}.[1-\frac{\left|\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|\right|}{\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|}\left]\right)---\left(2\right)$

定义4邻域形变主轴方向向量是用来描述像素点(p,q)邻域的局部各向异性形变中，最剧烈的形变方向和大小，记为

图像的邻域形变主轴方向向量则构成了该图像的邻域形变主轴向量场，即为Tagging MRI图像的邻域伸缩运动图像记为M_f。

根据上述定义，对心脏Tagging MRI图像的形变位移场M_u进行分析，获取图像的邻域形变主轴向量场也即图像的邻域伸缩运动图像，来估计纤维可能的走形方向。步骤如下：

Step1：根据定义1，求出心脏Tagging MRI图像中每个像素(p,q)的邻域相对形变场R_pq；若Tagging MRI图像所对应的时刻处于心脏收缩期，则转向step2；否则，若Tagging MRI图像所对应的时刻处于心脏舒张期，则转向step3。

Step2：根据每个像素(p,q)的邻域相对形变场R_pq，按照定义2(公式(1))计算该像素的邻域收缩主轴方向向量

获得图像的邻域收缩主轴方向场M_c；转Step4。

Step3：根据每个像素(p,q)的邻域相对形变场R_pq，按照定义3(公式(2))计算该像素的邻域拉伸主轴方向向量

获得图像的邻域收缩主轴方向场M_s；转Step4。

Step4：若心肌在该时刻处于收缩期内，则该Tagging MRI图像所对应的邻域形变主轴向量场M_f＝M_c；反之，若心肌在该时刻处于舒张期内，则M_f＝M_s。

通过(2)得到了Tagging MRI图像的邻域伸缩运动图像M_f，得到的邻域伸缩运动图像形如图4所示。

(3)计算(2)中获取的邻域伸缩运动图像M_f的运动中心点。运动中心点的计算过程如下：

过图像M_f上的每一像素点作一条垂直于该像素点形变方向的直线，计算运动中心点S₁到该直线的距离，同时将该距离乘以相应像素点形变大小的幅值，将所有像素点计算得到的该值相加，使得该值最大化的点o，即为所求的运动中心点，记图像M_f的运动中心点为S₁。如图5所示虚线是运动中心点到相应像素形变方向的距离，点o即为所求的运动中心点记为S₁。

(4)通过原始心脏DTI图像可以计算出每一个体素的扩散张量矩阵，张量矩阵表征了该体素中水分子扩散的各向异性方向和程度。由于水分子沿纤维束方向的扩散运动最快，所以一般认为水分子的主扩散方向平行于纤维束的走行方向。通过张量矩阵，可以计算出对应的主特征向量，就知道了该体素内水分子的主扩散方向，也即纤维束走向方向，通过这种方式获取了心脏DTI图像的心肌纤维走向图像M_a；

根据心肌纤维DTI图像计算获得纤维的走向图像，获取纤维走向图像的计算方法有流线跟踪法如Lori N F,Akbudak E,Shimony J S,et al.Diffusion tensor fiber tracking of human brain connectivity:aquisition methods,reliability analysis and biological results.NMRin Biomedicine,2002,15(7-8):494-515，张量弯曲法如Lazar M,Weinstein D M,Tsuruda J S,et al.White matter tractography usingdiffusion tensor deflection.Human brain mapping,2003,18(4):306-321等方法，记该纤维走向图像为M_a，纤维走向图像M_a形如图6所示。

(5)计算纤维走向图像M_a的中心点，记该中心点为S₂，该中心点的计算过程同(3)。

(6)通过平移运动将图像M_a与图像M_f的中心点S₂与运动中心点S₁对齐。

(7)运动中心点S₁与中心点S₂对齐以后，以中心点为圆心向外引若干条射线，这些射线将图像M_f与图像M_a等分成了大小相同的若干个扇形区域，记这些扇形的夹角大小为θ。设扇形的个数为n，将图像M_f的扇形区域分别记为x₁,x₂,x₃,…,x_n，将图像M_a的扇形区域分别记为y₁,y₂,y₃,…,y_n,图M_f与图M_a的扇形区域形如图7所示。

(8)从图像M_f与图像M_a的运动中心点S₁与S₂引若干条射线将每个扇形区域等分，将每个扇形区域内像素点形变方向投影到扇形平分线的垂线上，计算该扇形区域内所有投影值之和。设扇形区域x₁,x₂,x₃,…,x_n对应的投影和值分别为a₁,a₂,…,a_n，扇形区域y₁,y₂,y₃,…,y_n对应的投影和值分别为b₁,b₂,…,b_n。

(9)计算相似度度量函数，设相似度度量函数为f，f的计算过程如下：

设 $A=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{a}_i}{n},B=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_i}{n}$

相似度度量函数 $f=\frac{E\left\{\right[a_i-A]*[b_i-B\left]\right\}}{\sqrt{E\left\{{[a_i-A]}^2\right\}}*\sqrt{E\left\{{[b_i-b]}^2\right\}}}$

其中符号E{}表示求期望值。

(10)判断图像M_a旋转的角度是否穷举结束，即图像M_a从最初的位置每次顺时针方向旋转一个扇形角度的大小，判断图像M_a是否旋转了n次，n为图像M_a扇形的个数，旋转的过程如图7所示，内圆为M_a，外圆为M_f，图7是图像M_a旋转一个扇形角度过程的示意图，若图像M_a旋转的次数没有达到n次则转向步骤11，否则转向步骤12。

(11)将图像M_a顺时针旋转一个扇形大小的角度θ，并转向步骤(8)。当然此步骤中也可以沿逆时针方向旋转。

(12)通过将图像M_a的n次旋转，每次旋转后都能计算一个相似度度量函数f的值，找出使得f值最大时的旋转角度，并将该角度记为α。

(13)将图像M_a以中心点S₂为中心放大不同值的K倍，设依次放大K₁,K₂,…,K_m倍。每次放大后计算放大后图像与图像M_f重叠部分的相似度度量函数的值，设此时的相似度度量函数为g，其计算过程如下：

$g=\underset{i=1}{\overset{t}{\mathrm{\Σ}}}\left|\stackrel{\→}{u_i}\right|.\left|\stackrel{\→}{v_i}\right|.\cos (\stackrel{\→}{u_i},\stackrel{\→}{v_i})$

其中t为图像M_f与图像M_a重合的点的个数，

与

为两幅图像中相同位置的两像素点对应的方向向量，

与

为对应方向向量的模，

表示两向量夹角的余弦值。

通过m次放大能够得到m个不同的g值，找出使得相似度度量函数g最大化的K值，并记该K值对应的最大化相似度度量函数为g_t。

(14)至此配准一次结束，计算前后两次配准相似度度量函数最优值g_t与g_t-1之差的绝对值，判断该值是否小于给定的阈值r₁，同时判断扇形的角度是否小于给定的阈值r₂。若上述两个条件之一成立则配准搜索计算过程结束，同时转向步骤(16)，否则转向步骤(15)继续搜索。

(15)减小图像M_f与图像M_a对应的等分扇形角度的大小，并转向步骤(7)。

(16)通过上述各个步骤得到了图像M_a的平移参数，图像M_a的最佳旋转角度，图像M_a的最优放大系数，将图像M_a按照这些参数进行变换，即得到最后的配准图像。

本实施例中，是以纤维走向图像M_a相对邻域伸缩运动图像M_f进行平移、旋转及放大，当然也可以是邻域伸缩运动图像M_f相对纤维走向图像M_a进行平移、旋转及放大。

实施例：

本发明提出的心肌质点运动图像与心肌纤维走向图像配准方法涉及到若干参数，这些参数要针对具体处理的数据特点进行综合调节设定以达到最高配准精度，此处列出针对本发明处理数据集合而设定的参数：

步骤(7)将邻域伸缩运动图像M_f与纤维走向图像M_a以运动中心点为圆心等分为n个扇形，选取的n值为18。

步骤(13)在确定将图像M_a以中心点S₂为中心放大若干倍数时，均匀的在范围[0.5,2]内选择10个放大比例系数,即m的取值为10。

步骤(14)前后两次配准目标函数最优值g_t与g_t-1之差的绝对值，判断该值是否小于阈值r₁，或扇形的角度是否小于阈值r₂，阈值r₁取0.001，阈值r₂取4。

本发明方法是一种心脏多模态图像配准方法，能够解决心脏DTI图像与心脏Tagging MRI图像在纤维层次上的配准问题，DTI图像能够反映离体心脏的纤维结构，但不能反映活体心脏的运动信息，通过将心脏DTI图像与能提供活体心脏运动信息的Tagging MRI图像进行配准，来估计心肌纤维的运动信息。本发明方法根据心肌纤维具有沿着纤维伸展方向进行自主收缩与舒张运动的特点，将纤维的方向和心肌自主运动的方向作为物理含义截然不同而又密切相关的一对待配准特征，能够提供一种将DTI图像和Tagging MRI图像在纤维层次上配准的方法，为建立心脏纤维运动力学模型、在纤维层次上评估心脏运动功能打下基础，进而为分析心脏功能、心脏病病理学研究提供客观可靠的依据。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种心肌质点运动图像与心肌纤维走向图像配准方法，其特征在于，该方法包括下述步骤：

第1步，根据Tagging MRI图像序列获得心脏心肌质点的形变位移场图像M_u，也即为心肌的质点运动图像；

第2步，根据上述形变位移场图像M_u计算邻域伸缩运动图像M_f；

第3步，计算上述邻域伸缩运动图像M_f的运动中心点S₁；

第4步，获得原始心脏DTI图像的心肌纤维走向图像M_a；

第5步，计算上述纤维走向图像M_a的中心点S₂；

第6步，将上述邻域伸缩运动图像M_f的运动中心点S₁与上述纤维走向图像M_a的中心点S₂对齐，并确定平移量；

第7步，以对齐后的中心点为圆心向外引多条射线，将上述邻域伸缩运动图像M_f和纤维走向图像M_a分成大小相同的多个扇形；

第8步，计算邻域伸缩运动图像M_f和纤维走向图像M_a之间的相似度度量函数值f，

其中符号E{}表示求期望值，

a_i表示邻域伸缩运动图像M_f中第i个扇形区域内所有像素点形变方向投影到该扇形平分线的垂线上的投影值之和,b_i表示纤维走向图像M_a中第i个扇形区域内所有像素点形变方向投影到该扇形平分线的垂线上的投影值之和,n为邻域伸缩运动图像M_f和纤维走向图像M_a中扇形的个数；

第9步，将邻域伸缩运动图像M_f和纤维走向图像M_a相对旋转一个扇形角度大小，并按第8步的方法计算旋转后邻域伸缩运动图像M_f和纤维走向图像M_a之间的相似度度量函数值f，直到邻域伸缩运动图像M_f和纤维走向图像M_a之间相对旋转一周；

第10步，当邻域伸缩运动图像M_f和纤维走向图像M_a之间相对旋转一周得到每个旋转角度下的相似度度量函数值f时，确定上述相似度度量函数值f中的最大值，以及上述相似度度量函数值f的最大值所对应的旋转角度α；

第11步，相对邻域伸缩运动图像M_f和纤维走向图像M_a之间的初始位置，将邻域伸缩运动图像M_f和纤维走向图像M_a相对旋转上述角度α，并将旋转后的邻域伸缩运动图像M_f和纤维走向图像M_a中的任一图像相对另一图像以中心点为中心依次放大K₁,K₂,…,K_m倍，计算每次放大后邻域伸缩运动图像M_f和纤维走向图像M_a之间的相似度度量函数值g，得到m个相似度度量函数值，并确定这m个相似度度量函数值中的最大值g_t，以及最大值g_t所对应的放大系数；其中，m为放大的次数，相似度度量函数值g具体计算表达式为：

$g=\underset{i=1}{\overset{t}{\mathrm{\Σ}}}\left|\stackrel{\→}{u_i}\right|.\left|\stackrel{\→}{v_i}\right|.\cos (\stackrel{\→}{u_i},\stackrel{\→}{v_i})$

其中t为邻域伸缩运动图像M_f和纤维走向图像M_a重合的点的个数，

与

分别为两幅图像中相同位置的两像素点对应的方向向量，

与

分别为对应方向向量的模，

表示求两向量

与

的夹角的余弦值；

第12步，配准搜索计算本次迭代结束，判断本次最大相似度度量函数值g_t与上次最大相似度度量函数g_t-1之差的绝对值是否小于相似度阈值r₁，同时判断扇形的角度是否小于角度阈值r₂，若上述两个条件之一成立则配准搜索计算过程结束，同时转向第14步，否则转向第13步继续搜索；

第13步，减小邻域伸缩运动图像M_f和纤维走向图像M_a对应的等分扇形角度的大小，并转向第7步；

第14步，根据前述步骤中确定的图像间平移量，以及最佳旋转角度和最优放大系数，将邻域伸缩运动图像M_f或纤维走向图像M_a按照这些参数进行转换，得到配准结果。

2.如权利要求1所述的方法，其特征在于，所述第1步具体为：根据Tagging MRI图像序列计算获得心脏心肌质点形变位移场图像也即心脏心肌质点运动图像，心肌在拉伸与收缩的过程中，像素点(p,q)的形变用二维向量来表示，其中x，y分别为该像素点在水平与垂直方向上的位移量，该向量被定义为形变向量，心脏Tagging MRI图像中所有像素的形变向量构成了该图像的形变位移场，记整个心脏Tagging MRI图像的形变位移场图像为M_u。

3.如权利要求1或2所述的方法，其特征在于，所述第2步具体包括：

Step1：若Tagging MRI图像所对应的时刻处于心脏收缩期，则转向step2；否则，若Tagging MRI图像所对应的时刻处于心脏舒张期，则转向step3；

Step2：根据每个像素(p,q)的邻域相对形变场R_pq，按照公式

$\stackrel{\→}{L_c}=\underset{\stackrel{\→}{L}}{\arg \max }\left(\frac{1}{2}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\right|\stackrel{\→}{r_k}|.\cos {\mathrm{\θ}}_{(\stackrel{\→}{r_k},\stackrel{\→}{L})}.[1+\frac{\left|\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|\right|}{\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|}\left]\right)$

计算该像素的邻域收缩主轴方向向量

获得图像的邻域收缩主轴方向场M_c；转Step4；其中

定义为：心脏Tagging MRI图像中像素点(p,q)的邻域相对形变场中，将相对形变向量(k=1,2,3,…K)向过中心像素任意方向的直线L投影，方向向量记为

其中朝向中心像素的投影向量长度之和最大的直线方向，被称为该像素的邻域收缩主轴方向，其方向向量记为

K为邻域像素点个数；邻域收缩主轴方向场M_c为：所有像素的邻域收缩主轴方向向量组成了该图像的邻域收缩主轴方向场；

为向量

与向量

的夹角；

和分别为中心像素点到点A^(k)、B^(k)的矢量，点A^(k)、B^(k)分别为第k个邻域像素的相对形变向量

投影到过原点的直线L上的投影向量的起点与终点；||符号表示求模；

Step3：根据每个像素(p,q)的邻域相对形变场R_pq，按照公式

$\stackrel{\→}{L_s}=\underset{\stackrel{\→}{L}}{\arg \max }\left(\frac{1}{2}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}\right|\stackrel{\→}{r_k}|.\cos {\mathrm{\θ}}_{(\stackrel{\→}{r_k},\stackrel{\→}{L})}.[1-\frac{\left|\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|\right|}{\left|\stackrel{\→}{{\mathrm{OA}}^{\left(k\right)}}\right|-\left|\stackrel{\→}{{\mathrm{OB}}^{\left(k\right)}}\right|}\left]\right)$

计算该像素的邻域拉伸主轴方向向量获得图像的邻域拉伸主轴方向场M_s；转Step4；其中

定义为：心脏Tagging MRI图像中像素点(p,q)的邻域相对形变场中，将相对形变向量

(k=1,2,3,…K)向过中心像素任意方向的直线L投影，方向向量记为

其中远离中心像素的投影向量长度之和最大的直线方向，被称为该像素的邻域拉伸主轴方向，其方向向量记为

K为邻域像素点个数；邻域拉伸主轴方向场M_s为：所有像素的邻域拉伸主轴方向向量组成了该图像的邻域拉伸主轴方向场；

为向量

与向量

的夹角；

和

分别为中心像素点到点A^(k)、B^(k)的矢量，点A^(k)、B^(k)分别为第k个邻域像素的相对形变向量投影到过原点的直线L上的投影向量的起点与终点；||符号表示求模；

Step4：若心肌在该时刻处于收缩期内，则该Tagging MRI图像所对应的邻域伸缩运动图像M_f＝M_c；反之，若心肌在该时刻处于舒张期内，则M_f＝M_s。

其中邻域相对形变场R_pq的计算方式如下：将心脏Tagging MRI图像中像素点(p,q)的第k个邻域像素记为(p,q)^(k)，该邻域点相对于像素点(p,q)的相对形变向量为两像素点形变向量的差，记为

所有邻域像素的相对形变向量组成该像素的邻域相对形变场，记为R_pq。

4.如权利要求1至3任一项所述的方法，其特征在于，所述第3步具体包括：过邻域伸缩运动图像M_f上的每一像素点作一条垂直于该像素点形变方向的直线，计算运动中心点S₁到该直线的距离，同时将该距离乘以相应像素点形变大小的幅值，将所有像素点计算得到的该值相加，使得该值最大化的点即为所求的运动中心点，记邻域伸缩运动图像M_f的运动中心点为S₁。

Images