一种三维影像迭代配准方法与流程

1.本发明涉及图像处理领域，尤其涉及一种三维影像迭代配准方法。


背景技术：

2.随着生物医学工程和计算机技术的发展，医学影像学为临床诊断提供了多模态的医学图像。在当今临床邻域中，这些医学图像正扮演着越来越重要的角色，不仅应用于医疗诊断，还可用于手术计划的制定、放射治疗计划的制定、病理变化的跟踪和治疗效果评价等各个方面。
3.医学影像配准技术则是医学影像技术中的一个重要技术环节，其本质是通过寻找一种(或一系列)空间变换关系，使两幅图像的对应点到达空间位置或解剖结构上的完全一致，使得两个影像数据上的所有解剖点，或者至少是所有具有诊断意义的点都达到最佳匹配状态。
4.而现有技术中的配准方法存在配准误差较大、配准效率低等技术问题。


技术实现要素：

5.发明目的：本发明针对上述不足，提供了一种基于截面图像模板匹配的三维影像的迭代配准方法，可快速完成两个椎节的三维影像间的精确配准，效率高、适应性强，对精度较差的影像数据也能达到较好的配准效果。
6.技术方案：
7.一种三维影像迭代配准方法，包括步骤：
8.(1)获取基准影像和目标影像，并分别计算得到基准影像和目标影像中每一椎节的标准中心坐标，并据此得到基准影像和目标影像之间的初始变换关系；
9.(2)在基准影像中生成穿过各椎节中心的多个基准截面图像，根据步骤(1)得到的初始变换关系将基准影像和目标影像中对应的每一椎节的标准中心对准并分别采用前述多个基准截面图像在目标影像上进行截取得到对应的多个目标截面图像；
10.(3)将某一目标截面图像与对应的基准截面图像通过最佳相似匹配得到该目标截面图像的变换量，并据此计算得到该目标截面图像对应的目标截面坐标系在基准影像中的位姿；
11.(4)分别判断多个目标截面坐标系的位姿增量是否均小于设定阈值，若小于，则结合步骤(1)和步骤(3)计算得到基准影像坐标系和目标影像坐标系之间的最终变换关系；否则，分别更新多个目标截面图像的变换量及对应的目标截面坐标系的位姿，重复本步骤，直到多个目标截面坐标系的位姿增量均小于设定阈值，并据此得到基准影像坐标系和目标影像坐标系之间的最终变换关系。
12.所述步骤(1)中对应影像中每一椎节的标准中心坐标获取过程如下：
13.(11)计算得到基准影像和目标影像中每个椎节的初始中心坐标：
14.(12)将影像中的所有点对步骤(11)中计算出的每一椎节的初始中心坐标进行距
离聚类得到对应点集，并计算各点集的当前中心坐标；
15.(13)将影像中所有点对当前中心坐标进行距离聚类得到新的相应点集，并计算各新点集的当前中心坐标；
16.(14)重复步骤(13)，直至各点集的当前中心坐标不变为止，并将最终得到的当前中心坐标作为对应的各椎节的标准中心坐标。
17.所述步骤(11)具体如下：
18.已知影像在影像坐标系中的坐标范围为x
min
、x
max
、y
min
、y
max
、z
min
、z
max
，其中，椎节数量为n；
19.则得到每个椎节的初始中心坐标：
[0020][0021]
其中，p
i0
表示第i节椎节的初始中心点坐标，i表示影像中的第i节椎节。
[0022]
所述步骤(12)的距离聚类具体为：分别计算影像中每个点至每个椎节的初始中心坐标的距离，并将该点归类于距离其最近的初始中心坐标所对应的椎节，遍历影像中所有点，得到n个椎节点集。
[0023]
所述步骤(2)对基准影像和目标影像进行配准为：以步骤(1)得到的基准影像的标准中心坐标和目标影像的标准中心坐标为参考标记点，进行最小二乘刚体配准，得到基准影像和目标影像的变换关系。
[0024]
所述步骤(3)的最佳相似匹配具体如下：
[0025]
(31)在基准截面图像上设定基准roi图像，在目标影像的对应目标截面图像上遍历得到对应的目标roi图像并依次将其与基准roi图像进行相似性匹配评价，得到满足设定条件的第一匹配度量值对应的目标roi图像，并据此得到基准截面图像和对应目标截面图像的平移变换关系；
[0026]
(32)将目标截面图像对应的目标截面坐标系绕某一垂直于该截面图像坐标系且穿过该坐标系的原点的轴线以设定旋转角度为设定步长旋转，在每次旋转后计算目标截面图像与基准截面图像之间的匹配度量值，并在其小于第一匹配度量值时取前一次计算得到的匹配度量值作为第二匹配度量值；
[0027]
(33)改变旋转方向，以设定步长减小0.1
°
作为新步长将目标截面图像坐标系进行旋转，并在每次旋转后更新计算对应的目标截面图像与相应的基准截面图像之间的匹配度量值并比较其与第二匹配度量值的大小，直到其大于第二匹配度量值，此时得到目标截面图像的最佳旋转角度。
[0028]
第一匹配度量值为步骤(31)中计算得到的最大匹配度量值。
[0029]
第一匹配度量值为步骤(31)中计算得到的第一个大于0.9的匹配度量值。
[0030]
所述步骤(31)中，两幅截面图像的相似性匹配表示如下：
[0031]
设基准roi图像和目标roi图像分别为i1和i2，其中，基准roi图像定义为基准截面图像的中间区域图像，目标roi图像是通过在目标截面图像中设定目标roi图像框并从目标截面图像遍历截取得到；和分别表示i1、i2图像上的点p
ij
的灰度值，和分别表
示i1、i2图像上所有点的平均灰度值，则匹配度量值计算如下：
[0032][0033]
其中，η为匹配度量值。
[0034]
所述目标roi图像框遍历目标截面图像的步长设为1个像素值。
[0035]
所述步骤(3)具体为：
[0036]
针对某一椎节，设当前目标截面坐标系位姿下，最佳旋转角度为θ，则旋转变换为t
rot
；位移变换量分别为offx、offy、offz，则平移变换为t
trans
，则得到：
[0037][0038][0039]
所述步骤(4)具体为：
[0040]
所述多个目标截面图像分别为矢截面图像、轴截面图像和冠截面图像，设目标影像坐标系至基准影像坐标系的变换为t
transform
，其初始值为步骤(1)得到的初始变换关系；所述多个目标截面图像对应的矢截面坐标系、轴截面坐标系和冠截面坐标系的初始位姿分别为sa
tag0
、ss
tag0
、sc
tag0
，则当前的调整关系及各目标截面坐标系位姿调整为下式：
[0041][0042]
式中，分别表示当前目标截面图像的旋转变换和平移变换；sp
tag
依次取分别代入上式进行迭代计算，且在初次计算时代入初始位姿sa
tag0
、ss
tag0
、sc
tag0
进行计算；表示上一次迭代后目标影像至基准影像的变换关系，在初次计算时，代入初始变换值t0进行计算。
[0043]
所述步骤(4)中，增量阈值设为0.001，单个目标截面坐标系的位姿增量λ定义为位姿增量矩阵t
diff
的行列式值，描述如下：
[0044][0045]
其中，θ为最佳旋转角度；offx、offy、offz分别为目标roi图像匹配下位移变换量。
[0046]
有益效果：本发明可快速完成两个ct椎节影像间的精确配准，效率高，适应性强，对较差的影像数据也能达到较好的配准效果。
附图说明
[0047]
图1为本发明的流程图；
[0048]
图2为本发明得到的椎节中心示意图；
[0049]
图3为本发明中基准影像的截面图，其中，图3(a)为基准影像的矢截面示意图，图3(b)为基准影像的轴截面示意图，图3(c)为基准影像的冠截面示意图；
[0050]
图4为本发明中基准影像的截面示意图；
[0051]
图5为本发明中目标影像的截面图，其中，图5(a)为目标影像的矢截面示意图，图5(b)为目标影像的轴截面示意图，图5(c)为目标影像的冠截面示意图；
[0052]
图6为本发明的迭代配准流程图；
[0053]
图7为本发明中基准截面图像示意图；
[0054]
图8为本发明中基准截面图像与目标截面图像的相似性匹配的过程示意图。
具体实施方式
[0055]
下面结合附图和具体实施例，进一步阐明本发明。
[0056]
图1为本发明的流程图，如图1所示，本发明主要解决两个椎节ct影像的配准问题，两个椎节ct影像分别为基准影像和目标影像；本发明的方法包括以下步骤：
[0057]
(1)获取基准影像和目标影像，并分别建立对应的影像坐标系；
[0058]
(2)对影像中的椎节进行初步定位，以得到基准影像和目标影像中的每个椎节的标准中心坐标分别pi和pj，具体流程为：
[0059]
(21)计算基准影像和目标影像中每个椎节的初始中心坐标为p
i0
、p
j0
：设基准影像在基准影像坐标系中的坐标范围为x
min
、x
max
、y
min
、y
max
、z
min
、z
max
，坐标值范围变化最大的为脊椎的长度方向；每个影像中含有n节椎节，本发明以每个影像中含有三节椎节为例对具体步骤进行详细阐述，在其他实施例中，每个影像中椎节的数量根据实际情况而定；基准影像中每个椎节的初始中心坐标p
i0
采用下式进行计算：
[0060][0061]
其中，p
i0
表示第i节椎节的初始中心点坐标，i表示基准影像中的第i节椎节；
[0062]
采用同样的方法，计算出目标影像中每个椎节的初始中心坐标p
j0
，其中，j表示目标影像中的第j个椎节；
[0063]
(22)计算基准影像和目标影像中每个椎节的标准中心坐标pi、pj；
[0064]
(a)将基准影像中所有像素点对步骤(21)中计算出的对应椎节的初始中心坐标p
i0
进行聚类分类，即计算基准影像中每个像素点至每个椎节的初始中心坐标p
i0
的距离d
ik
＝|p
ik
p
i0
|，p
ik
表示基准影像中的第k个像素点；将该像素点p
ik
归类于距离其最近的初始中心坐标所对应的椎节，则可将基准影像中所有像素点分割到对应的三个椎节中的一个，得到三个椎节点集qi＝p
i1
,
…
,p
ik
,
…
,p
ik
；其中，k表示基准影像中像素点的数目；
[0065]
(b)通过三个椎节点集分别计算出每个椎节的当前中心坐标m为大于等于0的整数；
[0066]
(c)重复将基准影像中所有像素点对最新的当前中心坐标p
im
进行距离聚类并计算新的当前中心坐标p
im
，直至椎节的当前中心坐标不发生任何变化为止，并以最终的当前中心坐标p
im
作为每个椎节的标准中心坐标pi，如图2所示；
[0067]
(d)采用步骤22(a)、22(b)、22(c)的方法，计算出目标影像中每个椎节的标准中心坐标pj；
[0068]
(3)获取基准影像和目标影像的截面图像；
[0069]
(31)采用基准影像坐标系作为基准坐标系，并以基准影像中每个椎节的标准中心坐标作为原点建立每个椎节对应的截面图像的二维截面坐标系s
ai
、s
si
、s
ci
，分别称为第一矢截面坐标系、第一轴截面坐标系和第一冠截面坐标系，并在每个截面坐标系上分别截取尺寸为521mmx512mm的截面图像，形成基准截面图像，分别为第一矢截面图像、第一轴截面图像和第一冠截面图像，如图3(a)、3(b)、3(c)；其中，第一轴截面坐标系平行于基准影像坐标系的xy平面、第一矢截面平行于基准影像坐标系的yz平面，第一冠截面平行于基准影像坐标系的xz平面，基准影像的截面图像如图4所示；
[0070]
(32)以步骤(2)计算出的基准影像的标准中心坐标pi和目标影像的标准中心坐标pj为参考标记点，进行最小二乘刚体配准，完成基准影像和目标影像的整体初步配准得到二者的初始变换关系t0，并将此变换关系应用于目标影像，并采用步骤(31)中形成的第一矢截面图像、第一轴截面图像和第一冠截面图像分别截取目标影像形成对应的目标截面图像，分别为第二矢截面图像、第二轴截面图像和第二冠截面图像，如图5(a)、5(b)、5(c)，第二矢截面图像、第二轴截面图像和第二冠截面图像所在的三个二维坐标系分别为s
aj
、s
sj
、s
cj
，；
[0071]
(4)如图6所示，通过目标影像的截面图像与对应的基准影像的截面图像进行相似性匹配得到二者的变换关系，并更新当前目标截面坐标系位姿；
[0072]
具体为：
[0073]
(41)对目标截面图像与相应的基准截面图像进行相似性匹配得到最佳平移变换关系；
[0074]
基准截面图像与目标截面图像的相似性匹配的过程如图8所示，具体如下：
[0075]
(a)在基准截面图像上设定基准roi图像，基准roi图像定义为基准截面图像的中间区域图像且为正方形，基准roi图像宽度为基准截面图像宽度的1/2；在目标截面图像中设定目标roi图像框，目标roi图像框从目标截面图像的左上角开始进行遍历，依次截取目标roi图像框得到多个目标roi图像，直至目标roi图像框遍历至目标截面图像右下角；其中
遍历步长设为一个像素值；
[0076]
(b)依次将基准roi图像与步骤41中(a)步骤截取得到的目标roi图像进行相似性匹配评价计算，得到匹配度量值η和目标roi图像的位置偏移量；
[0077]
其中，两幅roi图像的相似性匹配可采取经典的相关系数评价值来评价，设基准roi图像和目标roi图像分别为i1和i2，和分别表示i1、i2图像上的点p
ij
的灰度值，和分别表示i1、i2图像上所有点的平均灰度值，则匹配度量值可表示为下式：
[0078][0079]
其中，η即为匹配度量值，η值越接近于1表示两幅截面图像越相似，越接近-1表示两幅截面图像越不相似参见图7；
[0080]
(c)选取步骤41中(b)步骤中匹配度量值最大时对应的基准roi图像和目标roi图像作为目标截面图像与基准截面图像的最佳相似性匹配结果，得到对应的平移匹配度量值为η1及对应的目标截面图像的偏移量；在其他实施方式中，为提高配准效率，可以在η大于等于0.9时停止相似性匹配评价计算，并以当前的基准roi图像和目标roi图像作为目标截面图像与基准截面图像的最佳相似性匹配结果；
[0081]
(d)根据步骤41中(c)步骤得到目标截面图像与相应的基准截面图像的平移匹配度量值η1及其对应的偏移量，即得到目标截面坐标系转换至基准截面坐标系下的偏移量offx，offy，offz，从而得到二者之间的平移变换t
trans
；
[0082]
(42)对目标截面图像与相应的基准截面图像进行相似性匹配迭代得到最佳旋转变换关系；
[0083]
(a)规定目标截面图像所对应的目标截面坐标系围绕某一轴顺时针或逆时针旋转的方向为角度增加方向，本发明中，该轴线垂直该截面图像坐标系且穿过该坐标系的原点，以设定的旋转角度作为设定步长，将该目标截面坐标系绕该轴线以设定的步长旋转，并更新目标截面图像；
[0084]
(b)每次旋转后，采用步骤(41)中(b)步骤的公式，更新计算该目标截面图像与相应的基准截面图像之间的第一旋转匹配度量值η2；比较当前的第一旋转匹配度量值η2与之前的平移匹配度量值η1的大小，若η2＞η1，则继续按当前旋转角度增加方向以设定步长继续绕该轴线旋转，直到η2＜η1，则取前一次计算得到的旋转匹配度量值为最终的第一旋转匹配度量值η2，转至步骤42中(c)步骤；
[0085]
(c)改变旋转方向，以设定步长减小0.1
°
作为新步长，将该目标截面坐标系绕该轴以新步长旋转，并在每次旋转后，更新目标截面图像，采用步骤(41)(b)步骤的公式，更新计算该目标截面图像与相应的基准截面图像之间的第二旋转匹配度量值η3并比较其与第一旋转匹配度量值η2的大小，直到η3＞η2，此时得到目标截面图像与相应的基准截面图像之间的最佳旋转匹配度量值η3，进而得到最佳旋转角度；
[0086]
(43)某一当前目标截面坐标系位姿更新按下式进行：
[0087]
针对某一椎节，设某一当前目标截面坐标系位姿下，当前目标截面图像的最佳旋转角度为θ，即旋转变换为t
rot
；目标roi图像匹配下位移变换量分别为offx、offy、offz，即
平移变换为t
trans
，则得到：
[0088][0089][0090]
根据上式对目标截面坐标系进行变换得到当前目标截面坐标系的位姿；
[0091]
(5)计算目标影像至基准影像的变换关系及第一矢截面坐标系、第一轴截面坐标系和第一冠截面坐标系的当前位姿：
[0092]
设目标影像至基准影像的变换关系为t
transform
，即为目标影像坐标系与基准影像坐标系的变换关系，其初始变换值为t0；第二矢截面坐标系、第二轴截面坐标系和第二冠截面坐标系的初始位姿分别为sa
tag0
、ss
tag0
、sc
tag0
，由步骤(31)中得到；则在当前调整下，当前目标影像至基准影像的变换关系及各目标截面坐标系位姿调整采用下式计算：
[0093][0094]
式中，分别表示当前目标截面图像的旋转变换和平移变换；sp
tag
依次取分别代入上式进行迭代计算，且在初次计算时代入初始位姿sa
tag0
、ss
tag0
、sc
tag0
进行计算；表示上一次迭代后目标影像至基准影像的变换关系，在初次计算时，代入初始变换值t0进行计算。
[0095]
(6)更新迭代目标影像坐标系的位姿，直到目标影像的三个目标截面坐标系的位姿增量小于设定阈值：
[0096]
针对单个椎节的三个目标截面坐标系的位姿采用步骤(4)进行更新之后，需依次计算目标影像中的三个目标截面坐标系(矢截面坐标系、轴截面坐标系、冠截面坐标系)的位姿增量；
[0097]
比较三个目标截面坐标系的位姿增量和设定阈值的大小，如果三个目标截面坐标系的位姿增量均小于阈值，则结束，并以当前的变换关系作为最终变换关系；如果三个目标截面坐标系的位姿增量中任一大于设定阈值，则采用步骤(5)的公式计算得到最新的三个目标截面坐标系的位姿并在目标影像中对三个目标截面坐标系的位姿进行更新，同时采用步骤(4)和步骤(5)的公式计算得到更新后的截面坐标系的平移变换t
trans
和旋转变
换t
rot
及变换关系据此计算三个目标截面坐标系的最新位姿增量并与设定阈值进行比较；重复本步骤，直到三个目标截面坐标系的位姿增量均小于设定阈值。
[0098]
本发明中，增量阈值设为0.001，单个截面坐标系的位姿增量λ定义为位姿增量矩阵t
diff
的行列式值，描述如下：
[0099][0100]
以上详细描述了本发明的优选实施方式，但是本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换(如椎节的数量等)，这些等同变换均属于本发明的保护范围。