技术领域
本发明涉及一种C型臂断层成像(C-armCT)系统中的离线标定方法，更具体地说，它 涉及一种基于螺旋线模体直接确定C-armCT系统中几何标定参数的方法，属于医疗成像领 域。
背景技术
集成了平板探测器的C型臂断层成像(C-armCT)系统，通过获得一系列二维X射线投 影数据去重建物体的断层图像。C型臂成像系统主要由X射线源，平板探测器，C型臂构架 组成。相对于传统CT而言，C-armCT具有移动方便，价格便宜等优点，尤其是其开放性 设计，能广泛地应用在介入治疗中的术前规划，术中检测及结果评估中。
但是在实际应用中，因为机械制造，系统集成及重力下垂等因素的影响，C-armCT系 统的实际轨道会偏离预设轨道。这些不规则的偏差会降低图像分辨率，或导致在重建图像中 出现严重伪影，影响医生判断。
目前针对此问题，螺旋线模体由于具有良好的3维空间信息，被广泛地应用在实际的几 何标定中。此模体主体为中空的圆柱形的树脂，上面嵌有按螺旋线均匀排列的一系列小钢珠， 这些小钢珠即为标志点。研究者采用的方法大多基于投影矩阵，即利用模体上标志点的3维 空间信息及投影图像的2维信息的对应关系来标定系统。但此方法没有对成像系统进行具体 地描述，如X射线源位置，探测器的倾斜及旋转角度，源到探测器的距离等，并且在系统偏 离较大的情况下，以此方法为基础得到重建图像质量可能会不理想。
另外一些研究者能获得系统的具体标定参数，但是此时螺旋线模体不再适用，需要设计 专门的标定模体，这样加大了额外工作量。
发明内容
本发明的目的是为了解决上述问题，提出一种基于螺旋线模体直接确定C型臂断层成像 系统中几何标定参数的方法。该方法利用常用的螺旋线模体，用解析方法获取了成像系统的 全部9个几何标定参数，具有较高的鲁棒性及准确性。
本发明利用螺旋线的几何特性，将螺旋线上的标志点按不同的分类方法分类，对每一组 而言，每四个标志点定义一个平行四边形。在不考虑投影点重叠的情况下，这四个标志点的 投影点确定一个四边形。这个投影得到的四边形对角线的交点即为平行四边形的对角线的交 点的投影。具体的技术方案如下：
第一步，将螺旋线模体放在C型臂置物台上，用C型臂扫描，获取投影数据。
第二步，确定螺旋线上标志点在平板探测器上的投影位置。
第三步，利用螺旋线的几何特性，确定模体坐标系z轴，x轴及两条平行于y轴的投影。
第四步，确定模体坐标系原点在平板探测器上的投影。
第五步，获取中间参数:a、b、c；
第六步，通过中间参数获得系统的9个几何标定参数：探测器的旋转角度（3个参数： 俯仰角，偏转角，倾斜角），X射线源到探测器距离（1个参数），X射线源位置（3个参数:x 方向，y方向，z方向），模体坐标系原点在探测器上的投影（2个参数：水平方向和竖直方 向）。
本发明的优点在于：
（1）该方法能获得系统的全部9个几何标定参数，从而能具体地描述此成像系统，对 机械制造有指导意义。并且获得独立的几何标定参数，为后期将其集成到相应的重建算法中 提供了可能性；
（2）该方法基于广泛应用的螺旋线模体，不需要额外加工标定模体；
（3）该方法是一种解析方法，不基于数学约化。并且充分利用螺旋线模体上标志点的信 息，保证了方法的鲁棒性及准确性。
附图说明
图1是螺旋线模体示意图。
图2是确定几何标定参数的流程图。
图3是确定坐标系z轴投影的示意图。
图中：
具体实施方式
下面将结合附图和实施例对本发明作进一步的详细说明。
本发明是一种确定C型臂断层成像系统中几何标定参数的方法，如图2所示，包括以下 几个步骤：
第一步，将螺旋线模体放在C型臂断层成像系统中的置物台上，用C型臂扫描，获取螺 旋线模体在平板探测器上的投影图像。
螺旋线模体如图1所示，本发明中模体镶嵌小钢珠(即标志点)的树脂部份跟通用的不大 一致，本发明采用螺旋管道代替了常用的中空圆柱形，这样就降低了制造模体所需要的树脂， 从而降低了制造此模体的成本（采用快速成型的加工方法）。小钢珠的排列跟通用的螺旋线状 模体是一致的，这样不会影响本方法的通用性。
螺旋线模体的放置要保证螺旋线的中心轴与C型臂系统旋转轴的方向大概一致。用C型 臂连续扫描一周，角度范围设定根据机型的不同，适用范围的不同调整。满足进行断层成像 要求的最小角度范围为180°＋最大扇角的一半。
最后，得到x射线穿透螺旋线模体在平板探测器上的投影图像。
第二步，确定螺旋线上标志点在平板探测器上的投影位置。
为了确定标志点投影在探测器上的位置，首先用canny边界检测的方法获得标志点的边 界，接着采用基于奇异值分解的椭圆拟合算法确定投影图像上标志点中心。本发明采用椭圆 的中心作为投影中心，根据相关文献，这两者间的误差可以忽略。
第三步，利用螺旋线的几何特性，确定模体坐标系z轴、x轴及两条平行于y轴的轴的 投影。
螺旋线上标志点的排列满足半径为R，螺距为2πh的螺旋线方程，第l个标志点位置的定 义如下：
rl＝(Rsinψl，Rcosψl，hψl)
ψl＝l×Δψ(l＝-N1，-N1+1，...，N2-1，N2)
其中：Δψ是螺旋线上标志点的角度间隔，由方程Δψ＝π/N给出。N的取值为整数，理 论上来说，N的取值范围可以从1到无穷大。但在实际应用中，当螺旋线长度保持固定，随 着N取值的增大，标志点的数量也随着增多，此时标志点的投影可能会在投影图像上出现重 叠，不利于后期处理。当N取值太小时，标志点数目过少，会影响算法精度。所以，在实际 应用中，N的取值范围为3≤N≤18是比较合理的。
参数N1和N2的取值为整数，它们的选取跟螺旋线中轴线的总长度和参考点的选取有关。 假设均匀排列的标志点构成的螺旋线中轴线的长度为L，那么N1+N2＝L/(hΔψ)-1。一般 选取位于螺旋线中部或接近中部的标志点做为参考点，此参数点的角度为0，即ψ0＝0。在给 出螺旋线中轴线的总长度及参考点的选取位置后，可以确定出N1和N2的取值。
（1）确定坐标系z轴投影的方法，如图3所示，从螺旋线上选取4个标志点，形成平行四边形如下：
P k z = { r A r B r C r D k | ψ A = ψ k , ψ B = ψ k + 2 π , ψ C = ψ k + 3 π , ψ D = ψ k + π , ]]>
k∈[-N1，N2-3N]}，
其中，rA，rB，rC，rD是螺旋线中的4个标志点，也是构成此平行四边形的4个顶点；ψA，ψB， ψC，ψD分别为点rA，rB，rC，rD在螺旋线中的角度位置；k是此平行四边形的索引值， ψk＝k×Δψ。
平行四边形对角线的交点为：
r I z k = ( 0,0 , h ( ψ k + 3 π / 2 ) ) . ]]>
平行四边形在探测器上的投影为四边形对角线交点在探测器上的投影为因为对角线交点在z轴上，所以投影得到的对角线交点在z轴的投影z′上，从新选取平行四边形，得到平行四边形对角线交点在探测器上的投影，根据两个投影对角线交点，确定模体坐标系z轴的投影z′轴。
图3是确定坐标系z轴投影的示意图，其中pA，pB，pC，pD分别为点rA，rB，rC，rD在探测器 上的投影点。在投影点不重叠的情况下，这4个投影点确定一个四边形，此四边形对角线的 交点为pI。pI与rI是一一对应关系。在图3中，设定螺旋线的角度间隔为π/6（，平行四边形 四个顶点rA，rB，rC，rD的在螺旋线中的角度参数分别为：
ψA＝-2π/3，ψB＝4ψ/3，ψC＝7ψ/3，ψD＝π/3，此时平行四边形对角线的交点坐标为对角线的交点在z轴上，因此投影得到的四边形的对角线交点确定了z轴的投影z′在探器上的一个点。同样的，可以得到一系列的投影四边形从而确定z轴的投影。
（2）确定坐标系x轴投影。从螺旋线上选取4个标志点，形成平行四边形如下：
P i x = { r A r B r C r D i | ψ A = ψ i , ψ B = ψ i + 2 π , ψ C = - ψ i , ψ D = - ψ i - 2 π , ]]>
i∈[max(-N1，-N2，-N+1)，min(N2-2N，N1-2N)]，i≠0}
其中，rA，rB，rC，rD是螺旋线中的4个标志点，也是构成此平行四边形的4个顶点；ψA，ψB，ψC，ψD分别为点rA，rB，rC，rD在螺旋线中的角度位置；i是此平行四边形的索引值，ψi＝i×Δψ。平行四边形对角线的交点为：
r I x i = ( R cos ψ i , 0,0 ) ]]>
与求z轴投影的类似，确定x轴上的一个点，得到相应的投影四边形对角线的交点，在选取平行四边形，得到平行四边形对角线交点在探测器上的投影，根据两个投影对角线交点，确定模体坐标系x轴的投影x′轴确定探测器上x轴的投影。
（3）确定两条平行于y轴的投影。从螺旋线上选取8个标志点，形成两个平行四边形和
P j 1 y = { r A r B r C r D j 1 | ψ A = ψ j 1 , ψ B = ψ j 1 + 2 π , ψ C = - ψ j 1 + π , π D = - ψ j 1 - π , ]]>
j 1 ∈ [ max ( - N 1 , N - N 2 , - N 2 + 1 ) , min ( N 1 - N , N 2 - 2 N ) ] , j 1 ≠ N 2 } , ]]>
P j 2 y = { r A r B r C r D j 2 | ψ A = ψ j 2 , ψ B = ψ j 2 + 2 π , ψ C = - ψ j 2 + π , π D = - ψ j 2 - π , ]]>
j 2 ∈ [ max ( 2 N - N 1 , N - N 2 ) , min ( N 1 - N , N 2 , N 2 - 1 ) ] , j 2 ≠ - N 2 } , ]]>
P j y = P j 1 y ∪ P j 2 y . ]]>
其中，rA，rB，rC，rD是螺旋线中的4个标志点，也是构成此平行四边形的4个顶点；ψA，ψB， ψC，ψD分别为点rA，rB，rC，rD在螺旋线中的角度位置；j1、j2是此平行四边形的索引值，
ψj1＝j1×Δψ,ψj2＝j2×Δψ。
两个平行四边形的交点为：
r I y j 1 = ( 0 , R sin ψ j 1 , πh / 2 ) ]]>
r I y j 2 = ( 0 , R sin ψ j 2 , - πh / 2 ) ]]>
对每个交点而言，它们的投影点用来确定两条平行于y轴的轴的投影。
与求z轴投影的类似，第一个平行y轴的y′轴投影：确定y轴上的一个点，得到相应的投影四边形对角线的交点，在选取平行四边形，得到平行四边形对角线交点在探测器上的投影，根据两个投影对角线交点，确定模体坐标系y轴的投影y′轴确定探测器上y轴的投影；
第二个平行y轴的y′′轴投影：确定y轴上的一个点，得到相应的投影四边形对角线的交点，在选取平行四边形，得到平行四边形对角线交点在探测器上的投影，根据两个投影对角线交点，确定模体坐标系y轴的投影y′′轴确定探测器上y轴的投影；
第四步，确定模体坐标系原点在平板探测器上的投影；
透入点为模体坐标系原点在平板探测器上的投影；模体坐标系原点在探测器上的投影由 x轴投影和z轴投影的交点确定。
对所有投影点的坐标进行预处理，将投影点的坐标减去透入点坐标。
第五步，获取中间参数a、,b、c；
ax,ay,az分别表示a在x、y、z轴的投影；bx,by,bz分别表示b在x、y、z轴的投影； cx,cy,cz分别表示c在x、y、z轴的投影；
对点其中：
u I z k = a z ( ψ k + 3 π / 2 ) h c z ( ψ k + 3 π / 2 ) h + 1 ]]>
v I z k = b z ( ψ k + 3 π / 2 ) h c z ( ψ k + 3 π / 2 ) h + 1 ]]>
其中：表示对应的在探测器上的坐标值。
是未知数为az，bz，cz的线性方程组。在k≥2的情况下，上述公式是有3个未知量，2k个 线性方程的超定线性方程组，能很容易通过例如基于奇异值分解的优化方法求解。
同样的，对于点其中：
u I x i = a x R cos ψ i c x R cos ψ i + 1 ]]>
v I x i = b x R cos ψ i c x R cos ψ i + 1 ]]>
其中，表示对应的在探测器上的坐标值。
ax，bx，cx可以通过上述方程组获得；
同样的，对于点其中：
u I y j 1 = a y R sin ψ j 1 + a z ( π / 2 ) h c y R sin ψ j 1 + c z ( π / 2 ) h + 1 ]]>
v I y j 1 = b y R sin ψ j 1 + b z ( π / 2 ) h c y R sin ψ j 1 + c z ( π / 2 ) h + 1 ]]>
u I y j 2 = a y R sin ψ j 2 + a z ( π / 2 ) h c y R sin ψ j 2 + c z ( π / 2 ) h + 1 ]]>
v I y j 2 = b y R sin ψ j 2 + b z ( π / 2 ) h c y R sin ψ j 2 + c z ( π / 2 ) h + 1 ]]>
其中：表示对应的在探测器上的坐标值；表示对应的在探测器上的坐标值；
ay，by，cy可以通过上述方程组获得。
第六步，通过中间参数获得C型臂断层成像系统的9个几何标定参数；
几何标定参数包括探测器的旋转角度、X射线源到探测器距离、X射线源的位置、原点 在探测器上的投影；
探测器的旋转角度包括俯仰角、偏转角、倾斜角，
X射线源的位置包括X射线源分别在x方向、y方向、z方向的位置；
模体坐标系原点在探测器上的投影包括水平方向投影和竖直方向投影。
通过上述得到的9个标定参数，对C型臂断层成像系统进行描述，对机械制造有指导作 用。大部份的C型臂重建算法根据投影矩阵来重建图像，在振动较大的时，此种方法不再适 用。获得独立的几何标定参数，能将此独立的几何参数集成到相应的重建算法中，获得较投 影矩阵更优的图像质量。
中间参数到9个几何标定参数的转换关系可以参考文献（C.Mennessier,R. Clackdoyle,andF.Noo,“Directdeterminationofgeometricalignmentparameters forcone-beamscanners,”PhysMedBiol,vol.54,no.6,pp.1633-60,Mar21, 2009.）
在对本发明的数值验证实验中，在添加相当于10%探测器像素大小的高斯噪声时，所测 得9个几何标定参数误差的均值：位置<0.02mm，角度<0.003°；误差的方差：位置<0.3mm, 角度<0.1°。可以看出本发明方法有较高的鲁棒性及准确性。