本发明公开了一种改进的CBCT几何参数标定方法，通过空间几何关系找到不同坐标系下的变换关系，以参数估计的方法提高标定参数的精确值。解析几何和特征点迭代优化相结合的方式精确计算标定参数，迭代估计方法的引入在提高精度的前提下也克服现有算法在标定过程对标定体模位置的要求，达到忽略体模位置而精确标定参数的目的。
1.一种改进CBCT几何参数标定方法，其特征在于，该方法的步骤如下：
1)按要求设计体模，体模中钢球位置误差在±0.001毫米以内；
2)360度全扇扫描，按顺时针或逆时针方向依次采集体模的投影图像；
3)对投影图像进行特征点检测与追踪，每一帧图像中检测的特征点形成上下两个椭圆，按顺时针方向对特征点坐标进行重新排序，再以差帧的方法追踪特征点，保证特征点一一对应，并在每副图像中利用检测到的特征点坐标拟合成两个椭圆，椭圆方程为：a(x-u0)2+b(y-v0)2+2c(x-u0)(y-v0)＝1；
其中，a、b、c为椭圆参数，u0、v0为椭圆中心坐标；
5)计算实际平板相对虚拟平板之间的旋转角度及初步计算
7)利用在步骤5)、6)中计算的参数关系，计算世界坐标系和投影图像之间的变换关系；
其中记 代表点Pp与点Pq的距离，A(p，q)为直线p，q之间夹角；P1，P2分别为两个椭圆中心；Pa是射线源在虚拟平面的投影点，其值可由以下公式获得：
其中ak、bk、ck、εk分别为特征点拟合的椭圆参数；其中Pφ＝(uφ,vφ)为两个椭圆长轴所在直线交点，与现实平板绕Y轴旋转角度有关；Pθ＝(uθ,vθ)为椭圆最两侧点连线构成的交点，与现实平板绕X轴旋转角度有关，其变换矩阵为：
5-3)计算射线源在各个坐标系中的坐标位置，通过以下计算公式可以粗略计算出射线源在各个坐标系中的位置：其中 分别代表射线源在虚拟坐标系和体模坐标系下的位置；
6-2)设以体模为基准，图像旋转角度为t，由投影关系可得射线源到体模坐标点的向量和射线源到投影点的坐标的向量对应成比例，得到两个比例关系；每一帧图像由特征点为已知量和步骤5)中得到的θ，φ，η为参数，使用Levenberg-Marquardt非线性优化算法进行参数估计， 以步骤5)中计算的结果为初值，t以上一帧计算的结果为初值，估计精确值；
步骤7)中计算世界坐标系和投影图像之间的变换关系的具体步骤如下：
7-1)每一帧计算一个射线源位置坐标，在步骤6)中计算的相对体模旋转的前提下，得到射线源在360度旋转扫描下的运动轨迹；
7-2)以第一帧射线源方向为参考方向，射线源轨迹将形成一个偏心圆，圆心即为等中心点；
7-3)利用已知参数，以等中心点为中心构成新的世界坐标系，进行坐标系参数变换。
一种改进CBCT几何参数标定方法
技术领域
[0001] 本发明涉及医学物理技术领域，具体涉及一种改进CBCT几何参数标定方法。
背景技术
[0002] CBCT图像引导系统主要用于辅助引导质子肿瘤治疗，通过与计划CT图像进行配准获取患者摆位误差。因此CT图像和患者的位置坐标需要精确定位。
[0003] 现有的标定算法大多假设标定体模旋转，而投影装置固定不变。但是对于CBCT装置此类方法在实际应用中存在先天上的缺陷。符合此类要求的CBCT标定算法，通过设计特定的标定体模，计算各种情况下参数变化，但是对于实际场景考虑不够全面，部分参数计算存在误差。
发明内容
[0004] 为了解决上述的技术问题，本发明的目的在于提供一种改进CBCT几何参数标定方法。
[0005] 本发明所要解决的技术问题为：
[0006] (1)如何进行部分参数的优化。
[0007] (2)如何适配实际应用场景。
[0008] 本发明的目的可以通过以下技术方案实现：一种改进CBCT几何参数标定方法，该方法的步骤如下：
[0009] 1)按要求设计体模，体模中钢球位置误差在±0.001毫米以内；
[0010] 2)360度全扇扫描，按顺时针或逆时针方向依次采集体模的投影图像；
[0011] 3)对投影图像进行特征点检测与追踪，每一帧图像中检测的特征点形成上下两个椭圆，按顺时针方向对特征点坐标进行重新排序，再以差帧的方法追踪特征点，保证特征点一一对应，并在每副图像中利用检测到的特征点坐标拟合成两个椭圆，椭圆方程为：
[0012] a(x-u0)2+b(y-v0)2+2c(x-u0)(y-v0)＝1；
[0013] 其中，a、b、c为椭圆参数，u0、v0为椭圆中心坐标；
[0014] 4)根据投影图像建立空间坐标系；
[0015] 5)计算实际平板相对虚拟平板之间的旋转角度及初步计算
[0016] 6)利用已知参数和初值进一步优化估计得到精确解；
[0017] 7)利用在步骤5)、6)中计算的参数关系，计算世界坐标系和投影图像之间的变换关系。
[0018] 进一步的，所述步骤5)的具体计算步骤如下：
[0019] 5-1)面内转角η利用如下公式计算：
[0020]
[0021] 其中记 代表点Pp与点Pq的距离，A(p，q)为直线p，q之间夹角；P1，P2分别为两个椭圆中心；Pa是射线源在虚拟平面的投影点，其值可由以下公式获得：
[0022]
[0023] 5-2)使用以下公式计算θ和φ角度：
[0024]
[0025]
[0026]
[0027]
[0028]
[0029] 其中ak、bk、ck、εk分别为特征点拟合的椭圆参数；其中Pφ＝(uφ,vφ)为两个椭圆长轴所在直线交点，与现实平板绕Y轴旋转角度有关；Pθ＝(uθ,vθ)为椭圆最两侧点连线构成的交点，与现实平板绕X轴旋转角度有关，其变换矩阵为：
[0030]
[0031] 5-3)计算射线源在各个坐标系中的坐标位置，通过以下计算公式可以粗略计算出射线源在各个坐标系中的位置：
[0032]
[0033]
[0034] 其中 分别代表射线源在虚拟坐标系和体模坐标系下的位置。
[0035] 进一步的，步骤6)中得到精确解的具体步骤如下：
[0036] 6-1)特征点位于虚拟坐标系下的坐标可表示为：
[0037]
[0038] 6-2)设以体模为基准，图像旋转角度为t，由投影关系可得射线源到体模坐标点的向量和射线源到投影点的坐标的向量对应成比例，得到两个比例关系；每一帧图像由特征点为已知量和5)中得到的θ，φ，η为参数，使用Levenberg-Marquardt非线性优化算法进行参数估计， 以步骤5)中计算的结果为初值，t以上一帧计算的结果为初值，估计精确值。
[0039] 进一步的，步骤7)中计算世界坐标系和投影图像之间的变换关系的具体步骤如下：
[0040] 7-1)每一帧计算一个射线源位置坐标，在步骤6)中计算的相对体模旋转的前提下，得到射线源在360度旋转扫描下的运动轨迹；
[0041] 7-2)以第一帧射线源方向为参考方向，射线源轨迹将形成一个偏心圆，圆心即为等中心点；
[0042] 7-3)利用已知参数，以等中心点为中心构成新的世界坐标系，进行坐标系参数变换。
[0043] 本发明的有益效果：
[0044] 本发明针对现有技术进行改进，提高参数计算的准确性。通过特征点的位置关系，采用优化算法的方式，迭代优化部分参数，克服在标定时由于体模摆放位置的不确定而带来的计算误差。
附图说明
[0045] 下面结合附图对本发明作进一步的说明。
[0046] 图1是设计标定体模的结构图；
[0047] 图2是步骤4)中建立的三个几何坐标系；
[0048] 图3是步骤4)中描述的实际平板相对虚拟平板的旋转关系；
[0049] 图4是标定体模所产生的投影图像其中的一帧；
[0050] 图5坐标系变化后射线源和投影点的坐标轨迹。
具体实施方式
[0051] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0052] 请参阅图1-5所示，本实施例提供了本发明的目的可以通过以下技术方案实现：一种改进CBCT几何参数标定方法，该方法的步骤如下：
[0053] 1)按要求设计体模，体模中钢球位置误差在±0.001毫米以内；首先设计如图1所示的标定体模，由钢球组成的两个圆形环组成，圆环半径为115毫米，封装在一个中空的丙烯酸圆柱体中。每个环包含12个钢球，每个钢球半径为2.5毫米，两环的间距为210毫米。利用激光校准，使体模尽可能置于等中心点，圆环中心连线与旋转轴尽可能平行。
[0054] 2)360度全扇扫描，按顺时针或逆时针方向依次采集体模的投影图像；对标定体模进行全扇扫描，按顺序采集720帧标定物体的投影图像，在每一帧图像的计算中，建立如图2所示坐标系，平板的偏转角度如图3所示。每一帧图像进行去噪、分割、阈值化等处理后得到类似于图4所示的图像。以圆检测的方法，检测到图像中的圆心作为特征点。
[0055] 3)对投影图像进行特征点检测与追踪，每一帧图像中检测的特征点形成上下两个椭圆，按顺时针方向对特征点坐标进行重新排序，再以差帧的方法追踪特征点，保证特征点一一对应，并在每副图像中利用检测到的特征点坐标拟合成两个椭圆，椭圆方程为：
[0056] a(x-u0)2+b(y-v0)2+2c(x-u0)(y-v0)＝1；
[0057] 其中，a、b、c为椭圆参数，u0、v0为椭圆中心坐标；
[0058] 4)根据投影图像建立空间坐标系；
[0059] 5)计算实际平板相对虚拟平板之间的旋转角度及初步计算 所述步骤5)的具体计算步骤如下：
[0060] 5-1)面内转角η利用如下公式计算：
[0061]
[0062] 其中记 代表点Pp与点Pq的距离，A(p，q)为直线p，q之间夹角；P1，P2分别为两个椭圆中心；Pa是射线源在虚拟平面的投影点，其值可由以下公式获得：
[0063]
[0064] 5-2)使用以下公式计算θ和φ角度：
[0065]
[0066]
[0067]
[0068]
[0069]
[0070] 其中ak、bk、ck、εk分别为特征点拟合的椭圆参数；其中Pφ＝(uφ,vφ)为两个椭圆长轴所在直线交点，与现实平板绕Y轴旋转角度有关；Pθ＝(uθ,vθ)为椭圆最两侧点连线构成的交点，与现实平板绕X轴旋转角度有关，其变换矩阵为：
[0071]
[0072] 5-3)计算射线源在各个坐标系中的坐标位置，通过以下计算公式可以粗略计算出射线源在各个坐标系中的位置：
[0073]
[0074]
[0075] 其中 分别代表射线源在虚拟坐标系和体模坐标系下的位置。
[0076] 6)利用已知参数和初值进一步优化估计得到精确解；步骤6)中得到精确解的具体步骤如下：
[0077] 6-1)特征点位于虚拟坐标系下的坐标可表示为：
[0078]
[0079] 6-2)设以体模为基准，图像旋转角度为t，由投影关系可得射线源到体模坐标点的向量和射线源到投影点的坐标的向量对应成比例，得到两个比例关系；每一帧图像由特征点为已知量和5)中得到的θ，φ，η为参数，使用Levenberg-Marquardt非线性优化算法进行参数估计， 以步骤5)中计算的结果为初值，t以上一帧计算的结果为初值，估计精确值。
[0080] 7)利用在步骤5)、6)中计算的参数关系，计算世界坐标系和投影图像之间的变换关系。上文计算出所有坐标系之间的变换关系，然后将体模坐标系下的参数转化到以等中心点为中心的世界坐标系下的参数，得到转换坐标系后射线源和等中心点投影点在旋转360度下的运动轨迹如图5所示。
[0081] 步骤7)中计算世界坐标系和投影图像之间的变换关系的具体步骤如下：
[0082] 7-1)每一帧计算一个射线源位置坐标，在步骤6)中计算的相对体模旋转的前提下，得到射线源在360度旋转扫描下的运动轨迹；
[0083] 7-2)以第一帧射线源方向为参考方向，射线源轨迹将形成一个偏心圆，圆心即为等中心点；
[0084] 7-3)利用已知参数，以等中心点为中心构成新的世界坐标系，进行坐标系参数变换。
[0085] 以上内容仅仅是对本发明结构所作的举例和说明，所属本技术领域的技术人员对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离发明的结构或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。
