一种倾斜图像的迭代校正方法与流程

1.本发明涉及图像处理技术领域，尤其涉及一种倾斜图像的迭代校正方法。


背景技术：

2.锥形束ct技术是采用x射线光源和平板探测器对物体进行旋转扫描，获取物体在不同角度下的投影，然后运用锥形束反投影重建算法获取物体三维图像的技术。实现锥形束反投影重建算法，需要为其提供由x射线源、探测器和物体组成系统的几何参数。这些几何参数对重建图像的分辨率和伪影等有重要影响，为获取高质量的重建ct图像，需要精确标定这些几何参数。
3.现有几何参数标定方法主要依赖于二维图像信息，多采用分步求解的几何标定方法，这些标定方法往往没有考虑各个几何参数间的相互影响问题，如探测器在不同方向倾斜情况下的影响及ct设备在不同旋转角度间形变造成的相互影响；另外，现有标定方法多采用“汇聚点”计算倾角，但由于不同方向角度倾斜的影响，“汇聚点”的求解精度波动较大，从而影响整体的几何参数标定精度。


技术实现要素：

4.发明目的：为解决上述问题，本发明提供了一种倾斜图像的迭代校正方法，在保证精度的同时提高了角度求解的鲁棒性。
5.技术方案：
6.一种倾斜图像的迭代校正方法，包括步骤：
7.(1)通过ct对标定模体进行扫描得到标定模体在平板探测器上的投影图像，并获取图像中各钢珠的中心坐标，其中，钢珠设置为两层并平行排布在标定模体内，每层钢珠组成一个标准的圆形；
8.(2)分别计算得到投影图像上形成的椭圆的方程，并据此求解对应的椭圆参数；
9.(3)根据标定模体的设计参数及步骤(2)得到的椭圆参数及射线源至平板探测器距离计算得到当前平板探测器相对于其理想位置的当前倾角，其中射线源至平板探测器距离的初始值为设计值；
10.(4)以平板探测器当前倾角对图像进行倾角校正，并计算校正后椭圆参数及射线源到平板探测器的实际距离；
11.(5)根据标定模体的设计参数及步骤(4)得到椭圆参数及射线源到平板探测器的实际距离计算得到平板探测器校正后相对于其理想位置的校正后的倾角；
12.(6)判断步骤(5)中校正后的倾角是否小于预设阈值，若是，则结束循环；否则，以步骤(5)得到的校正后的倾角与步骤(4)中的当前倾角之和作为最新的当前倾角，返回步骤(4)。
13.所述步骤(1)中，所述标定模体中，每层钢珠的数量设置为4～16个，钢珠的直径设为1～5mm。
14.所述标定模体中，每个钢珠的直径相同。
15.每层钢珠的数量设置为12个，每个钢珠的直径设为2mm。
16.所述标定模体为圆柱体，两层钢珠对称排布在标定模体的中心的上下两侧。
17.所述步骤(1)中，获取其中各钢珠的中心坐标具体为：
18.对图像进行高斯滤波及阈值分割，并以钢珠投影的像素点个数确定分割阈值，通过自适应阈值分割算法初步提取钢珠的位置坐标，并据此对钢珠进行边缘提取分割，获取各钢珠的精确轮廓，采用最小二乘法进行椭圆拟合获取得到各钢珠中心坐标。
19.所述步骤(3)中计算平板探测器相对于其理想位置的倾角具体为：
20.平板探测器相对于其理想位置绕平行于其竖直边的轴线的倾角为φ，绕平行于其水平边的轴线的倾角为θ；
21.根据式(1)、(2)分别计算得到对应倾角φ
i
和θ
i
；
[0022][0023][0024]
其中，v
0φ
、v
1φ
为平板探测器相对于其理想位置绕平行于其竖直边的轴线的倾角为φ
i
的情况下，上层钢珠形成的圆在实际平板探测器上投影形成的椭圆的长轴的顶点；v
2φ
、v
3φ
为平板探测器相对于其理想位置绕平行于其竖直边的轴线的倾角为φ
i
的情况下，下层钢珠形成的圆在实际探测器上投影形成的椭圆的长轴的顶点；v
0θ
、v
1θ
为平板探测器相对于其理想位置绕平行于其水平边的轴线的倾角为θ
i
的情况下，上层钢珠形成的圆在实际探测器上投影形成的椭圆的长轴的顶点；v
2θ
、v
3θ
为平板探测器相对于其理想位置绕平行于其水平边的轴线的倾角为θ
i
的情况下，下层钢珠形成的圆在实际探测器上投影形成的椭圆的长轴的顶点；d2为射线源到平板探测器的距离，其初始值为设计值。
[0025]
所述步骤(3)中计算平板探测器相对于其理想位置的倾角还包括：以对应倾角φ
i
和θ
i
对投影图像进行倾斜校正，获取校正后图像内的椭圆参数，并据此计算平板探测器绕射线源在平板探测器投影与射线源连线为轴线的倾角η
i
。
[0026]
所述倾角η
i
的计算具体为：
[0027]
v
0φ
v
1φ
和v
2φ
v
3φ
所在直线的交点为p
φ
，p
φ
在平板探测器上的垂点为p
a
，由式(3)计算得到点p
a
至原点o的距离d
oa
：
[0028][0029]
其中，d
o1
为上层椭圆中心至原点o的距离，d
o2
为下层椭圆中心至原点o的距离，a1、b1为上层椭圆方程的系数，a2、b2为下层椭圆方程的系数；
[0030]
根据式(4)计算得到平板探测器绕射线源在平板探测器投影与射线源连线为轴线的倾角η
i
：
[0031]
η＝[d
1a
a(x，l1)+d
2a
a(x，l2)]/d
12
ꢀꢀꢀ
(4)
[0032]
其中，l1、l2分别表示v
0φ
v
1φ
和v
2φ
v
3φ
所在直线，a(x，l1)和a(x，l2)分别为和所在直线与平板探测器上平行于其水平边的轴线的夹角；d
1a
、d
2a
分别表示上下层椭圆中心至点pa的距离，d
12
表示上下层椭圆中心之间的距离。
[0033]
所述步骤(5)中计算校正后椭圆对应参数及射线源到实际平板探测器的距离具体为：
[0034]
由式(5)计算得到射线源s在当前平板探测器所在的实际坐标系中的坐标值：
[0035][0036]
根据下式计算得到射线源s到实际平板探测器的实际距离d2：
[0037][0038]
所述步骤(6)具体为：
[0039]
根据步骤(4)中的当前倾角计算此时射线源s到实际平板探测器的实际距离d2，并根据步骤(4)得到的椭圆参数结合式(1)、(2)(3)、(4)计算得到校正后的倾角φ
i
′
、θ
i
′
、η
i
′
；
[0040]
若φ
i
′
的绝对值小于第一阈值，则令φ
i
′
＝0；若φ
i
′
的绝对值不小于第一阈值，则令φ
i+1
＝φ
i
+φ
i
′
作为最新当前倾角，并对图像进行倾角φ
i
′
的校正；
[0041]
若θ
i
′
的绝对值小于第一阈值，则令θ
i
′
＝0；若θ
i
′
的绝对值不小于第一阈值，则令θ
i+1
＝θ
i
+θ
i
′
作为最新当前倾角，并对图像进行倾角θ
i
′
的校正；
[0042]
若η
i
′
的绝对值小于第二阈值，则令η
i
′
＝0；若η
i
′
的绝对值不小于第二阈值，则令η
i+1
＝η
i
+η
i
′
作为最新当前倾角，并对图像进行旋转角η
i
′
的校正。
[0043]
所述第一阈值设为δ＝0.00087rad，所述第二阈值设为δ1＝0.000087rad。
[0044]
步骤(5)中，校正后的倾角小于预设阈值时，根据当前倾角计算标定参数，并输出。
[0045]
有益效果：本发明依据透视模型理论，利用探测器在不同方向旋转引起物体投影大小的变化计算角度误差，在保证精度的同时提高了角度求解的鲁棒性。同时，为消除不同参数间的相互影响问题，本发明采用迭代标定，逐次校正的方式，在标定的过程中逐步校正几何参数误差，直至满足标定精度要求。
附图说明
[0046]
图1为x射线源、标定模体及平板探测器的投影示意图；
[0047]
图2为本发明的流程图；
[0048]
图3为平板探测器存在倾角φ时标定模体投影示意图；
[0049]
图4为平板探测器存在倾角θ时标定模体投影示意图；
[0050]
图5为平板探测器同时存在倾角φ、θ时标定模体投影示意图；
[0051]
图6为平板探测器存在倾角η时标定模体投影示意图。
具体实施方式
[0052]
下面结合附图和具体实施例，进一步阐明本发明。
[0053]
本发明揭示的实施例中，标定模体的设计参数如下：标定模体为圆柱形，在标定模体内平行排布两层钢珠，每层钢珠的数量为12，钢珠直径为2mm，每层12个钢珠安置在标准圆形轨道上，标准圆形轨道对称排布在标定模体的中心的上下两侧，使得两层钢珠对称排布在标定模体的中心的上下两侧，该标准圆形轨道的半径为r，单位mm；上下两个标准圆形
轨道之间的竖直距离为d，单位mm。本发明采用锥形束ct的x射线源扫描标定模体并将其投影至平板探测器，x射线源穿过标定模体，x射线源发出的射线的中心线垂直于平板探测器，但是实际操作过程中，平板探测器相较于其理想位置会有角度偏差，本发明通过迭代计算得到平板探测器相较于理想位置的旋转角度，以校正几何参数误差。
[0054]
具体的，在进行校正之前，将标定模体放置在视野中央，通过正位和侧位拍摄图像将其调整到视野的正中央，然后通过锥形束ct对标定模体进行旋转扫描，锥形束ct的旋转轴线与标定模体的中心线基本重合，得到序列图像，如图1所示，x射线源s发出x射线束穿过标定模体，投影在平板探测器上产生投影图像。上下两层钢珠形成的标准圆在平面探测器上的投影为椭圆，上下两个椭圆各自相距最远的顶点，即椭圆的长轴顶点分别为v0、v1、v2、v3。
[0055]
图2为本发明的倾斜图像的迭代校正方法的流程图，如图2所示，本发明的倾斜图像的迭代校正方法包括如下步骤：
[0056]
(1)通过锥形束ct对标定模体进行旋转扫描，获取得到标定模体在平板探测器上的投影图像；
[0057]
(2)获取各钢珠的中心坐标：对单张投影图像进行高斯滤波及阈值分割，并以钢珠投影的像素点个数确定分割阈值，通过自适应阈值分割算法初步提取钢珠的位置坐标，并据此对钢珠进行边缘提取分割，获取各钢珠的精确轮廓，采用最小二乘法进行椭圆拟合获取得到各钢珠中心坐标；
[0058]
(3)对步骤(2)得到的上下层钢珠按逆时针进行排序，并对上下层钢珠进行中心对称配对，则各对中心对称钢珠连线之间的交点即为标定体模的中心，并据此求解得到标定模体的中心点投影在平板探测器上的坐标；
[0059]
(4)对上下层钢珠集合分别进行最小二乘椭圆拟合，分别获得上下层钢珠的椭圆方程，求解得到上下层钢珠在平板探测器上投影组成的椭圆的长轴顶点及椭圆中心坐标；
[0060]
(5)计算平板探测器倾角；
[0061]
图3为探测器存在倾角φ时标定模体在平板探测器上的投影示意图，如图3所示，虚线矩形表示平板探测器的理想位置，实线矩形表示平板探测器的实际位置；在平板探测器处于理想位置时，x射线源s穿过标定模体的中心的射线在平板探测器上的投影为o；以o点为原点，以穿过o点并平行于平板探测器的竖直边的轴线为y轴，穿过原点o并平行于平板探测器水平边的轴线x轴，穿过o点且同时垂直于x轴和y轴的轴为z轴建立理想坐标系；
[0062]
则平板探测器相对于平板探测器的理想位置绕y轴的倾角为φ时，平板探测器所在的实际坐标系与理想坐标系之间存在夹角，此时实际坐标系由x
φ
轴、y
φ
轴及z
φ
轴组成，其中，x
φ
轴为x轴绕y轴对应方向旋转角度φ得到，y
φ
轴与y轴一致，z
φ
轴为z轴绕y轴对应方向旋转角度φ得到；
[0063]
图4为平板探测器存在倾角θ时标定模体在平板探测器上的投影示意图，如图4所示，虚线矩形表示理想情况下的平板探测器的理想位置，实线矩形表示平板探测器的实际位置；平板探测器相对于平板探测器的理想位置绕x轴的倾角θ时，平板探测器实际所在的实际坐标系与平板探测器的理想位置所在的理想坐标系之间存在夹角，此时实际坐标系由x
θ
轴、y
θ
轴与z
θ
轴组成，其中，x
θ
轴与x轴一致，y
θ
轴为y轴绕x轴旋转角度θ得到，z
θ
轴为z轴绕x轴旋转角度θ得到；
[0064]
图5为探测器同时存在倾角φ、θ时标定模体投影示意图，根据式(1)、(2)分别计算得到平板探测器存在的倾角φ和θ的当前倾角值φ
i
和θ
i
；
[0065][0066][0067]
其中，v
0φ
、v
1φ
为平板探测器相对于平板探测器的理想位置绕y轴的倾角为φ
i
的情况下，上层钢珠形成的圆在实际平板探测器上投影形成的椭圆的长轴的顶点；v
2φ
、v
3φ
为平板探测器相对于平板探测器的理想位置绕y轴的倾角为φ
i
的情况下，下层钢珠形成的圆在实际平板探测器上投影形成的椭圆的长轴的顶点；v
0θ
、v
1θ
为平板探测器相对于平板探测器的理想位置绕x轴的倾角为θ
i
的情况下，上层钢珠形成的圆在实际平板探测器上投影形成的椭圆的长轴的顶点；v
2θ
、v
3θ
为平板探测器相对于平板探测器的理想位置绕x轴的倾角为θ
i
的情况下，下层钢珠形成的圆在实际平板探测器上投影形成的椭圆的长轴的顶点；d1为x射线源s到标定模体中心的距离，d2为射线源s到平板探测器的距离，初次计算时，d1、d2为设计值，在后期迭代过程中，d1、d2的值为式(7)中计算得到的实际距离；
[0068]
以前述计算得到的当前平板探测器存在倾角φ
i
和θ
i
对投影图像进行倾斜校正，获取校正后图像内的椭圆参数，并据此计算平板探测器相对于平板探测器的理想位置绕z轴的倾角η的当前值η
i
，图6为平板探测器存在倾角为η时标定模体投影示意图，具体计算如下：
[0069]
如图3所示，v
0φ
v
1φ
和v
2φ
v
3φ
所在直线的交点为p
φ
，p
φ
在y
φ
轴上的垂点为p
a
，由式(3)可计算得到：
[0070][0071]
其中，d
oa
为点p
a
至原点o的距离，d
o1
为上层椭圆中心至原点o的距离，d
o2
为下层椭圆中心至原点o的距离，a1、b1为上层椭圆方程的系数，a2、b2为下层椭圆方程的系数；
[0072]
根据式(4)计算平板探测器相对于平板探测器的理想位置绕z轴的倾角η的当前值η
i
；
[0073]
η＝[d
1a
a(x，l1)+d
2a
a(x，l2)]/d
12
ꢀꢀꢀ
(4)
[0074]
其中，l1、l2分别表示v
0φ
v
1φ
和v
2φ
v
3φ
所在的直线，a(x，l1)和a(x，l2)分别为和所在直线与x轴的夹角；d
1a
、d
2a
分别表示上下层椭圆中心至点p
a
的距离，d
12
表示上下层椭圆中心之间的距离；
[0075]
(6)计算平板探测器被校正后，x射线源s到标定模体中心的实际距离d1和x射线源s到平板探测器的实际距离d2；
[0076]
根据平板探测器相对于平板探测器的理想位置存在倾角的当前倾角值φ
i
、θ
i
和η
i
，计算旋转中心q(即标定模体的中心)、x射线源s在世界坐标系和平板探测器实际所在的实际坐标系中的坐标值，并据此计算得到x射线源s到标定模体中心的实际距离d1和x射线源s到平板探测器的实际距离d2；具体为：
[0077]
(61)由式(5)计算得到x射线源s在当前平板探测器所在的实际坐标系中的坐标值：
[0078][0079]
同时由式(6)计算得到旋转中心q在当前平板探测器所在的实际坐标系中的坐标值：
[0080][0081]
(62)根据式(7)计算得到x射线源s到标定模体中心的实际距离d1和x射线源s到平板探测器的实际距离d2：
[0082][0083]
(7)计算经当前值φ
i
、θ
i
、η
i
校正后的平板探测器与其理想位置之间的校正后的倾角φ
i
′
、θ
i
、η
i
′
，并与设定阈值进行比较；
[0084]
根据式(7)计算得到的新的x射线源s到标定模体中心的实际距离d1和x射线源s到平板探测器的实际距离d2，结合式(1)、(2)分别计算得到校正后的平板探测器的校正后的倾角φ
i
′
、θ
i
′
；
[0085]
设定第一阈值δ＝0.00087rad，若φ
i
′
的绝对值小于第一阈值，则令φ
i
′
＝0；若φ
i
′
的绝对值不小于第一阈值，则令φ
i+1
＝φ
i
+φ
i
′
，并对图像进行倾角φ
i
′
进行校正；
[0086]
若θ
i
′
的绝对值小于第一阈值，则令θ
i
′
＝0；若θ
i
′
的绝对值不小于第一阈值，则令θ
i+1
＝θ
i
+θ
i
′
，并对图像倾角θ
i
′
进行校正；
[0087]
通过式(4)计算得到校正后的平板探测器的倾角η
i
′
，设定第二阈值δ1＝0.000087rad，若η
i
′
的绝对值小于第二阈值，则令η
i
′
＝0；若η
i
′
的绝对值不小于第二阈值，则令η
i+1
＝η
i
+η
i
′
，并对图像倾角η
i
′
进行校正；
[0088]
(8)判断φ
i
′
、θ
i
′
、η
i
′
是否都为0，若是，则结束循环，以当前倾角作为实际平板探测器与理想探测器位置之间的夹角计算标定参数并输出标定参数，标定参数为根据当前倾角计算得到的旋转中心、射线源在平板探测器坐标系中的坐标值、射线源到标定模体中心的实际距离d1和射线源到平板探测器的实际距离d2；否则，以φ
i+1
、θ
i+1
、η
i+1
作为最新当前值φ
i
、θ
i
、η
i
返回步骤(6)。
[0089]
在其他实施例中，每层钢珠的数量可以设置为4～16个，钢珠的直径可以设为1～5mm，且各钢珠中心在同一个圆周上；优选的，所有钢珠的直径相等。
[0090]
在本发明中，第一阈值和第二阈值的具体数值可以根据实际情况设定。
[0091]
本发明依据透视模型理论，利用平板探测器在不同方向旋转引起物体投影大小变化计算角度误差，在保证精度的同时提高了角度求解的鲁棒性。同时，为消除不同参数间的相互影像问题，本发明采用迭代标定、逐次校正的方式，逐步校正几何参数误差，直至满足标定精度要求。
[0092]
以上详细描述了本发明的优选实施方式，但是本发明并不限于上述实施方式中的
具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换(如数量、形状、位置等)，这些等同变换均属于本发明的保护范围。