CN111973204B - 一种纳入重力的新型双平板x光机的校准方法 - Google Patents
一种纳入重力的新型双平板x光机的校准方法 Download PDFInfo
- Publication number
- CN111973204B CN111973204B CN202010775024.XA CN202010775024A CN111973204B CN 111973204 B CN111973204 B CN 111973204B CN 202010775024 A CN202010775024 A CN 202010775024A CN 111973204 B CN111973204 B CN 111973204B
- Authority
- CN
- China
- Prior art keywords
- ray
- dimensional point
- lateral
- ray system
- flat
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 53
- 230000005484 gravity Effects 0.000 title claims abstract description 43
- 230000000694 effects Effects 0.000 title description 2
- 229910000831 Steel Inorganic materials 0.000 claims abstract description 61
- 239000010959 steel Substances 0.000 claims abstract description 61
- 238000012795 verification Methods 0.000 claims abstract description 22
- 238000005457 optimization Methods 0.000 claims abstract description 14
- 230000009466 transformation Effects 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 33
- 239000013598 vector Substances 0.000 claims description 16
- 238000010586 diagram Methods 0.000 claims description 12
- 238000013507 mapping Methods 0.000 claims description 10
- 238000003384 imaging method Methods 0.000 claims description 9
- 238000002360 preparation method Methods 0.000 claims description 9
- 238000000513 principal component analysis Methods 0.000 claims description 9
- 238000005516 engineering process Methods 0.000 claims description 8
- 238000001514 detection method Methods 0.000 claims description 7
- 230000009977 dual effect Effects 0.000 claims description 5
- 238000004458 analytical method Methods 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 230000003287 optical effect Effects 0.000 claims description 3
- 238000004064 recycling Methods 0.000 claims description 3
- 238000012800 visualization Methods 0.000 claims description 3
- 238000011160 research Methods 0.000 abstract description 5
- 230000005021 gait Effects 0.000 abstract description 2
- 239000002184 metal Substances 0.000 description 6
- 210000000988 bone and bone Anatomy 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 239000011521 glass Substances 0.000 description 3
- 230000000399 orthopedic effect Effects 0.000 description 3
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 241000287196 Asthenes Species 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003339 best practice Methods 0.000 description 1
- 235000019219 chocolate Nutrition 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000037230 mobility Effects 0.000 description 1
- 230000002980 postoperative effect Effects 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/58—Testing, adjusting or calibrating thereof
- A61B6/582—Calibration
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/58—Testing, adjusting or calibrating thereof
- A61B6/582—Calibration
- A61B6/585—Calibration of detector units
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- Engineering & Computer Science (AREA)
- Radiology & Medical Imaging (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- High Energy & Nuclear Physics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Length-Measuring Devices Using Wave Or Particle Radiation (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
本发明提供了一种纳入重力的新型双平板X光机的校准方法，使用相机模型来模拟X光***，采用直接线性变换求取参数的解析解，并将此解作为后续最小二乘问题的初解，使用列文伯格‑马夸尔特最优化算法继续迭代的方法，完成了双平板X光***的校准。为了方便后续基于双平板X光的步态研究，本发明引入了铅重力线以测量重力方向。为了验证校准精度，本发明使用单元块搭建了一个嵌有钢珠的验证盒，通过比较由校准参数还原出的钢珠三维坐标和厂家实测的三维坐标的误差，我们达到了均方根误差(root mean squared error，RMSE＝0.1847mm)的重建精度。
Description
技术领域
本发明涉及双平板X光机的校准技术，具体涉及一种纳入重力的新型双平板X光机的校准方法。
背景技术
双平面X光是骨科研究中的一种无创的高精度的成像方式。它在计算机辅助的术前规划，计算机智能诊断，解剖三维重建，生物力学运动追踪等领域都有重要应用。例如，骨科研究中，经常需要分析骨的运动学数据，进而比较正常骨与异常骨的运动模式，从而揭示骨科疾病的机理或者评估术后效果。双平面X光因为其无创、高精度等特点，已逐渐成为主流的运动追踪方式。双平面X光***包含两套成一定角度交叉放置的X光设备，拍摄时将感兴趣的骨置于两套X光锥形辐照区域重叠处，以此来模拟人的双目视觉从而进行骨的定位与追踪。
为了还原X光***的成像环境，需要通过特殊的校准手段(calibration)来测取：1.每一套X光自身的参数；2.两套X光间的位姿关系。掌握这些成像参数是所有基于双平面X光的后续应用的重要前置基础，因此其校准精度非常重要。
现有的X光成像设备主要有：1.基于影增(intensifier)的C臂机(c-armfluoroscopy)2.基于平板的X光机(planer digital imaging system)。这两种设备因为成像原理不同，设备大小不同，因此在组装成双平面X光时，所需要的配套校准技术也不同。由于C臂机X光在医院场景下更容易获取，且移动性强，因此组装成双平面X光***也更为容易，对应的校准技术也较为成熟。
基于C臂机的双平面X光的校准思路如下：制作一个参数已知的校准盒。校准盒有双层，每一层是一块塑料圆盘，圆盘内嵌入x光无法透过的钢珠，两层圆盘间用若干塑料支柱连接。校准时，将校准盒贴在影增屏上固定住，通过实际获取的图像结合已知的校准盒的物理信息，可以推算出单X光机的光源相对于影增的位置信息，也可以推算影增带来的各种图像扭曲的参数，从而完成对单平面的校准。
通过对两台C臂机分别执行以上操作，可分别对两台C臂机进行单独的校准。完成C臂机各自的校准后，需要确定两台C臂机的位置关系。方法是制作一透明玻璃板，玻璃板的四个角嵌有螺丝。让此玻璃板分别被两台C臂机照射到。由于每一台C臂机都校准过了，所以可以推算出物体和单C臂机的位置关系。由于是拍摄的同一物体，借助物体这一中间媒介，利用坐标转换的矩阵乘法，可以推算出双C臂机的位置关系。
C臂机双平面X光***由于自身物理硬件原因，导致对应的校准方法存在以下问题：
1.为了完成校准需要执行三步操作(两台X光机的分别校准，两台X光机之间的校准)，效率低。若两台C臂机相对位置关系移动，则需要重复校准。
2.校准盒以及校准方法是基于影增的C臂机开发设计的，但是对于平板X光机并不适用，因为平板相对影增要大很多，设计一个能够贴住平板的校准盒不太现实。
3.由于双C臂机每次的摆位都不同，所以重力相对于X光***的坐标表达每一次都会变，需要每次花大量时间去检测重力方向，进一步降低效率。
目前已有厂商开发出了基于平板X光机的双平面X光***，但是对应的校准方法仍旧是一片空白。本专利提供一种纳入重力方向的且适配双平板X光校准技术，希望为此领域提供最佳实践的范式。
发明内容
针对现有技术中的缺陷，本发明的目的是提供一种纳入重力的新型双平板X光机的校准方法。本发明的技术方案如下：
一种纳入重力的新型双平板X光机的校准方法，包括如下步骤：
S1：建立一套双平板X光机***，包括：
正位X光***、侧位X光***；
所述侧位X光***包括：第一平板X光机；
所述正位X光***包括：第二平板X光机；
所述第一平板X光机包括：第一X光源、第一平板探测器；
所述第二平板X光机包括：第二X光源、第二平板探测器；
所述第一平板X光机、第二平板X光机可固定地整体做上下移动；
S2：准备一校验盒；
所述校验盒包括依次相连的第一侧面、第二侧面、第三侧面、第四侧面；第一侧面和第三侧面对应第一平板X光机；第二侧面和第四侧面对应第二平板X光机；
所述第一侧面靠近第一平板探测器，第三侧面靠近第一X光源；将第一平板X光机拍摄的图片称为侧位图；
所述第二侧面靠近第二平板探测器，第四侧面靠近第二X光源；将第二平板X光机拍摄的图片称为正位图；
所述校验盒的每个侧面上均设置有若干个钢珠；且每个侧面的钢珠的分布方式一致；每个侧面的各个钢珠的三维坐标均已知；
S3：获取预存的校验盒的每个侧面各个钢珠的三维坐标；
S4:将校准盒摆在两台平板X光机光路重叠的区域，两台平板X光机同时曝光，获取正位图及侧位图；所述正位图和侧位图上均出现若干个黑点，每一个黑点均为一个对应钢珠的投影；
S5：进行侧位X光***和正位X***校准前的对应关系准备；
S6：求出侧位X光***的内、外参数以及正位X光***的内、外参数；
S7：利用公式(1)求取正位X光***坐标系、侧位X光***坐标系的位姿关系；如下：
其中，记Rl,tl为侧位X光***的外参数，Rf,tf为正位X光***的外参数；
tl＝-RlCl；
tf＝-RfCf；
Pcf,Pcl是同一个点分别在正位X光***坐标系和侧位X光***坐标系下的坐标；
S8:将侧视图和正视图分别利用计算机图形学进行展示，即使用计算机图形学分别将侧位X光***、正位X光***的虚拟环境重现出来；；
S9：确定重力方向，即确定重力向量在双平板X光***下的坐标表达。
可选地，所述步骤S5进一步包括：
S51：首先进行侧位X光***的校准准备工作，包括：
S511：使用“Hough圆检测”算法检测侧位图上黑点的圆心坐标，得到这些黑点的二维坐标，并将这些二维坐标集合在一起，记为二维点集xy，这里记录共有n个二维坐标；n为正整数；
S512：读取预存的钢珠坐标中涉及到侧位图的那一部分钢珠的坐标，将这些钢珠的坐标记为三维点集XYZ，这里记录共有n个三维坐标；
S513：建立从三维点集XYZ到二维点集xy的映射，得到二者的对应关系，共有n个对应关系。
S52：使用S51相同的方法，进行正位X光***的校准准备工作。
可选地，所述步骤S6进一步包括：
S61：首先求出侧位X光***的内、外参数，包括：
S611：建立空间三维点到图像二维点之间的变换矩阵：
从空间中三维点到图像上的二维点，借助齐次坐标，可用一次矩阵乘法完成：
其中：
Pw是步骤S51中三维点集XYZ的齐次形式，为4×1向量；p是像素坐标的齐次形式，为3×1向量；M是3×4矩阵，M矩阵蕴含了内参数和外参数；
S612：将第i个三维点集XYZ和二维点集xy的对应关系代入上式，并改写为齐次线性方程组的形式：
其中：
pi＝[ui vi]T是二维点集xy的第i个二维点；
Pi＝[Xi Yi Zi]T是三维点集XYZ中第i个三维点；
i＝1,…,n；
S613:将所有的n组对应关系联立，有如下过定线性齐次方程组：
记为
Am＝0
S614:对A奇异值分解：
A＝UDVT
U、D、V是A的三个分解矩阵；
过定方程的最优解是V矩阵的最后一列，即最优解
S615:构建非线性最优化问题：
其中：
pi＝[ui vi]T是二维点集xy的第i个二维点
Pi＝[Xi Yi Zi]T是三维点集XYZ中的第i个三维点坐标；
使用
作为此优化问题的初解，使用levernberg-marquardt最优化算法迭代求解，记优化终解为
S616:将向量
重排为矩阵形式，得到
S617:由
矩阵解析K,R,C；
其中：K是内矩阵，R是联系校准盒坐标系和侧位X光***坐标系的旋转矩阵，C是第一光源坐标；将K称为内参数，将R、C称为外参数；
将
矩阵左边3×3部分记为
则
对
使用QR分解，并将结果求逆可得K和R；
将
矩阵右边3×1部分记为
则
S62：使用S61相同的方法，求取正位X光***的内、外参数。
可选地，所述步骤S8进一步包括：
S81：将侧视图利用计算机图形学进行展示，即使用计算机图形学将侧位X光***的虚拟环境重现出来；其进一步包括：
S811：将全局坐标系选在侧位X光***的原点；
S812：根据侧位X光***的内参数虚拟出一个虚拟成像平面，然后用纹理贴图的办法将实际获取的侧位图贴图到这个虚拟平面上；
S813：根据据侧位X光***的外参数，将校验盒上涉及侧位成像的钢珠在空间中用黑色点显示出来；
S82：使用S81相同的方法，将正视图利用计算机图形学进行展示。
可选地，所述步骤S9进一步包括：
选取一重力方向的指征物，其同时被正位X光***、侧位X光***拍摄，分别获得所述指征物的正位图和侧位图；利用所述指征物在正位图和侧位图上的投影、以及之前求出的正位X光***、侧位X光***各自的内、外参数，推算重力方向。
可选地，所述步骤S9进一步包括：
选取一重力方向的指征物，并预先获得所述指征物的三维模型，使用正位X光***、侧位X光***对所述指征物进行拍摄，获得正位图、侧位图；
将预先获得的指征物的三维模型导入到正位X光***、侧位X光***的虚拟环境中；
将视点放在第二光源、第一光源处，同时手动调整所述指征物的三维模型在虚拟环境中的位置和朝向，当从视点处发现指征物的模型覆盖住了正位图、侧位图上指征物的轮廓，认为此时三维模型的位置和朝向信息还原出了指征物的实际位置和朝向信息；
根据指征物自身的形状特点、以及之前手动调整获得的指征物的实际空间位置和朝向信息，进一步推算重力方向信息。可选地，所述方法还包括：
步骤S10：准备一验证盒，在S2的校验盒拍摄完毕后，保持第一平板X光机、第二平板X光机的位置不变，用验证盒替换校验盒，并拍摄正位图和侧位图；结合步骤S6获得的正、侧位X光***的内、外参数，使用“三角测量”技术，验证重建三维点的精度。
可选地，所述步骤S10进一步包括：
S101：使用乐高积木搭成一25cm*25cm*25cm的正方体盒子作为验证盒，将27颗钢珠嵌入正方体盒子内；
所述钢珠以近似3*3*3的晶格的形式分布于正方体内；钢珠之间在两个方向上间隔相同，在第三个方向间隔不同；
所述验证盒还内嵌入用钢丝制成的十字，圆圈，方形和三角标志；正、侧位图上的黑点即对应上述钢珠；
S102:预存所有27颗钢珠的三维坐标点，作为预存三维点集；
S103：使用“Hough”圆检测算法检测正、侧位图上的黑点；
将侧位片上的点集记为S100,将正位片上的点集记为S200；
S104:使用“三角测量”技术：于S100中的每一点，向第一光源引出一条射线，再取S200中的每一点，向第二光源引出一条射线，计算两射线的空间距离；
选择能产生最小射线空间距离的点对，认为这一对点是同一空间三维点在两侧图像上的投影；
同时计算两射线最靠近时，最小二乘意义上的交点，作为重建点三维坐标；
循环使用此方法，可建立S100与S200的对应关系，对每一点对，也同时重建了其对应的三维点坐标；利用以上方法计算出来的27颗钢珠三维点坐标的集合，即为重建的三维点集；
S105:步骤S104的重建的三维点集是位于双平板X光***坐标系下的，步骤S102的预存三维点集是位于验证盒坐标系下的，故将步骤S102的预存点集和步骤S104的重建的三维点集配准到同一坐标系下，方法如下：
结合步骤S6获得的正、侧位X光***的内、外参数，对重建三维点集使用PCA确定其中心，从而确定两个坐标系的偏移向量；
PCA给出的前三个主成分分量其实就是旋转过的三坐标轴；
利用钢珠排列在两个方向上间隔相同在第三个方向间隔不同这一空间异质性，配合可视化来看是否需要坐标轴的翻转来更好的拟合两组点；
当从可视化中看出重建的三维点集和预存的三维点集基本重合后，再用暴力求解法确定点的对应关系；
最终可计算出RMSE来表征校准精度。
可选地，第一平板X光机、第二平板X光机的光轴存在一夹角，所述夹角为30度-150度。
与现有技术相比，本发明具有如下的有益效果：
1、本发明的校准只需拍摄一次即可，相比基于C臂机的双平面X光***的三次拍摄更加方便高效。
2、通过引入铅重力线，计算得到了重力方向在X光平板坐标系下的表达，为后续运动追踪和相关生物力学研究提供参考基础。
3、双平板***是锁死在移动杆上的，相对位置关系固定，重力方向相对于坐标系固定，所以***校准和重力方向的检测都可以只做一次就完成。
4、通过本方法可以精确在计算机中还原出实际的物理投影模型，是后续手动配准及自动计算的先导条件。验证的亚毫米级精度完全满足临床研究中的追踪精度要求，且为后续的三维重建和配准提供了精度保障。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1是本发明具体实施例一种纳入重力的新型双平板X光机的校准方法的流程图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
如图1，一种纳入重力的新型双平板X光机的校准方法，包括如下步骤：
S1：建立一套双平板X光机***，包括：
正位X光***、侧位X光***；
所述侧位X光***包括：第一平板X光机；
所述正位X光***包括：第二平板X光机；
第一平板X光机、第二平板X光机的光轴呈一夹角；所述夹角为30度-150度。夹角的范围这里仅为举例，本发明不对此做出限定。
所述第一平板X光机包括：第一X光源、第一平板探测器；
所述第二平板X光机包括：第二X光源、第二平板探测器；
所述第一平板X光机、第二平板X光机可相对固定地整体做上下移动；
S2：准备一校验盒；
所述校验盒包括依次相连的第一侧面、第二侧面、第三侧面、第四侧面；第一侧面和第三侧面对应第一平板X光机；第二侧面和第四侧面对应第二平板X光机；
所述第一侧面靠近第一平板探测器，第三侧面靠近第一X光源；将第一平板X光机拍摄的图片称为侧位图；
所述第二侧面靠近第二平板探测器，第四侧面靠近第二X光源；将第二平板X光机拍摄的图片称为正位图；
所述校验盒的每个侧面上均设置有若干个钢珠；且每个侧面的钢珠的分布方式一致；每个侧面的各个钢珠的三维坐标均已知；
S3：获取预存的校验盒的每个侧面各个钢珠的三维坐标；
S4:将校准盒摆在两台平板X光机光路重叠的区域，两台平板X光机同时曝光，获取正位图及侧位图；所述正位图和侧位图上均出现若干个黑点，每一个黑点均为一个对应钢珠的投影；
S5：进行侧位X光***和正位X***校准前的对应关系准备；
S6：求出侧位X光***的内、外参数以及正位X光***的内、外参数；
S7：利用公式(1)求取正位X光***坐标系、侧位X光***坐标系的位姿关系；如下：
其中，记Rl,tl为侧位X光***的外参数，Rf,tf为正位X光***的外参数；
tl＝-RlCl；
tf＝-RfCf；
Pcf,Pcl是同一个点分别在正位X光***坐标系和侧位X光***坐标系下的坐标；
S8:将侧视图和正视图分别利用计算机图形学进行展示，即使用计算机图形学分别将侧位X光***、正位X光***的虚拟环境重现出来；；
S9：确定重力方向，即确定重力向量在双平板X光***下的坐标表达。
S10：准备一验证盒，在S2的校验盒拍摄完毕后，保持第一平板X光机、第二平板X光机的位置不变，用验证盒替换校验盒，并拍摄正位图和侧位图；结合步骤S6获得的正、侧位X光***的内、外参数，使用“三角测量”技术，验证重建三维点的精度。
其中，所述步骤S5进一步包括：
S51：首先进行侧位X光***的校准准备工作，包括：
S511：使用“Hough圆检测”算法检测侧位图上黑点的圆心坐标，得到这些黑点的二维坐标，并将这些二维坐标集合在一起，记为二维点集xy，这里记录共有n个二维坐标；n为正整数；
由于钢珠排列规则，钢珠投影点在图像上以规则的栅格点方式排列，因此靠近校准盒侧面中心的钢珠自然对应图像中心的黑点，而远离侧面中心的钢珠对应图像边角的黑点。利用这种规则的排列规律和投影常识，可以建立从三位点集XYZ到二维点集xy的映射(即知道哪个三维点投影到哪个二维点)。
S512：读取预存的钢珠坐标中涉及到侧位图的那一部分钢珠的坐标，将这些钢珠的坐标记为三维点集XYZ，这里记录共有n个三维坐标；
S513：建立从三维点集XYZ到二维点集xy的映射，得到二者的对应关系，共有n个对应关系。
S52：使用S51相同的方法，进行正位X光***的校准准备工作。
其中，所述步骤S6进一步包括：
S61：首先求出侧位X光***的内、外参数，包括：
S611：建立空间三维点到图像二维点之间的变换矩阵：
从空间中三维点到图像上的二维点，借助齐次坐标，可用一次矩阵乘法完成：
其中：
Pw是步骤S51中三维点集XYZ的齐次形式，为4×1向量；p是像素坐标的齐次形式，为3×1向量；M是3×4矩阵，M矩阵蕴含了内参数和外参数；
S612：将第i个三维点集XYZ和二维点集xy的对应关系代入上式，并改写为齐次线性方程组的形式：
其中：
pi＝[ui vi]T是二维点集xy的第i个二维点；
Pi＝[Xi Yi Zi]T是三维点集XYZ中第i个三维点；
i＝1,…,n；
S613:将所有的n组对应关系联立，有如下过定线性齐次方程组：
记为
Am＝0
S614:对A奇异值分解：
A＝UDVT
U、D、V是A的三个分解矩阵；
过定方程的最优解是V矩阵的最后一列，即最优解
S615:构建非线性最优化问题：
其中：
Pi＝[Xi Yi Zi]T是三维点集XYZ中的第i个三维点坐标；
使用
作为此优化问题的初解，使用levernberg-marquardt最优化算法迭代求解，记优化终解为
S616:将向量
重排为矩阵形式，得到
S617:由
矩阵解析K,R,C；
其中：K是内矩阵，R是联系校准盒坐标系和侧位X光***坐标系的旋转矩阵，C是第一光源坐标；将K称为内参数，将R、C称为外参数；
将
矩阵左边3×3部分记为
则
对
使用QR分解，并将结果求逆可得K和R；
将
矩阵右边3×1部分记为
则
S62：使用S61相同的方法，求取正位X光***的内、外参数。
其中，所述步骤S8进一步包括：
S81：将侧视图利用计算机图形学进行展示，即使用计算机图形学将侧位X光***的虚拟环境重现出来；其进一步包括：
S811：将全局坐标系选在侧位X光***的原点；
S812：根据侧位X光***的内参数虚拟出一个虚拟成像平面，然后用纹理贴图的办法将实际获取的侧位图贴图到这个虚拟平面上；
S813：根据据侧位X光***的外参数，将校验盒上涉及侧位成像的钢珠在空间中用黑色点显示出来；
S82：使用S81相同的方法，将正视图利用计算机图形学进行展示。
其中，所述步骤S9进一步包括：
选取一重力方向的指征物，其同时被正位X光***、侧位X光***拍摄，分别获得所述指征物的正位图和侧位图；利用所述指征物在正位图和侧位图上的投影、以及之前求出的正位X光***、侧位X光***各自的内、外参数，推算重力方向。推算方法具体为：
在双平板X光机***拍摄校准图像时，于校准盒中心悬挂一铅垂线，用铅垂线方向来代表重力方向；(此处，铅垂线即为指征物。)
铅垂线于正位图和侧位图上各自留下了直线的投影；
利用正位图上的直线的投影、第二X光源，确定一个第一平面；
利用侧位图上的直线的投影、第一X光源，确定一个第二平面；
第一平面和第二平面在空间的交线，即为铅垂线的方向，也即重力方向。
需要说明的是，本实施例，确定重力方向的方法同时依靠正位X光***、侧位X光***两个***。具体实施时，仅选取一个X光***，也可以实现重力方向的确定，方法为：选取一重力方向的指征物，并预先获得所述指征物的三维模型，使用正位X光***、侧位X光***对所述指征物进行拍摄，获得正位图、侧位图；
将预先获得的指征物的三维模型导入到正位X光***、侧位X光***的虚拟环境(由步骤S8获得)中；
将视点放在第二光源、第一光源处，同时手动调整所述指征物的三维模型在虚拟环境中的位置和朝向，当从视点处发现指征物的模型覆盖住了正位图、侧位图上指征物的轮廓，认为此时三维模型的位置和朝向信息还原出了指征物的实际位置和朝向信息；
根据指征物自身的形状特点、以及之前手动调整获得的指征物的实际空间位置和朝向信息，进一步推算重力方向信息。具体的推算方法为：
定制一巧克力棒状的金属物作为指征物，用铅垂线悬吊起来，被正位X光***、侧位X光***拍摄。预先获取所述金属物的激光扫描模型(即三维模型)，并将金属物的三维模型导入正位X光***、侧位X光***的虚拟环境中手动调整位置和朝向。
利用所述金属物可看作一长方体形状，选取其某一边拟合直线，并利用前面获得的金属物的最终位置和朝向信息，将直线方向向量认为是重力方向。
其中，所述步骤S10进一步包括：
S101：使用单元块搭成一25cm*25cm*25cm的正方体盒子作为验证盒，将27颗钢珠嵌入正方体盒子内，实际操作时，选取对应的乐高积木，将钢珠嵌入乐高积木的槽内；本实施例中，单元块可以使用乐高积木。这里仅为举例，本发明不对此做出限定。
所述钢珠以近似3*3*3的晶格的形式分布于正方体内；钢珠之间在两个方向上间隔相同，在第三个方向间隔不同；
所述验证盒还内嵌入用钢丝制成的十字，圆圈，方形和三角标志；正、侧位图上的黑点即对应上述钢珠；
S102:预存所有27颗钢珠的三维坐标点，作为预存三维点集；
S103：使用“Hough”圆检测算法检测正、侧位图上的黑点；
将侧位片上的点集记为S100,将正位片上的点集记为S200；
S104:使用“三角测量”技术：于S100中的每一点，向第一光源引出一条射线，再取S200中的每一点，向第二光源引出一条射线，计算两射线的空间距离；
选择能产生最小射线空间距离的点对，认为这一对点是同一空间三维点在两侧图像上的投影；
同时计算两射线最靠近时，最小二乘意义上的交点，作为重建点三维坐标；
循环使用此方法，可建立S100与S200的对应关系，对每一点对，也同时重建了其对应的三维点坐标；利用以上方法计算出来的27颗钢珠三维点坐标的集合，即为重建的三维点集；
S105:步骤S104的重建的三维点集是位于双平板X光***坐标系(侧位X光***坐标系)下的，步骤S102的预存三维点集是位于验证盒坐标系下的，故将步骤S102的预存点集和步骤S104的重建的三维点集配准到同一坐标系下，方法如下：
结合步骤S6获得的正、侧位X光***的内、外参数，对重建三维点集使用PCA(Principalcomponentanalysis)确定其中心，从而确定两个坐标系的偏移向量；
PCA给出的前三个主成分分量其实就是旋转过的三坐标轴；
利用钢珠排列在两个方向上间隔相同在第三个方向间隔不同这一空间异质性，配合可视化来看是否需要坐标轴的翻转来更好的拟合两组点；
当从可视化中看出重建的三维点集和预存的三维点集基本重合后，再用暴力求解法确定点的对应关系；
最终可计算出RMSE(均方根误差，root mean squared error)来表征校准精度。
综上所述，本实施例使用相机模型来模拟X光***，采用直接线性变换(directlinear transform,DLT)求取参数的解析解，并将此解作为后续最小二乘问题的初解，使用列文伯格-马夸尔特最优化算法(levernberg-marquardt)继续迭代的方法，完成了双平板X光***的校准。为了方便后续基于双平板X光的步态研究，我们引入了铅重力线以测量重力方向。为了验证校准精度，我们使用乐高搭建了一个嵌有钢珠的小盒子，通过比较由校准参数还原出的钢珠三维坐标和厂家实测的三维坐标的误差，我们达到了均方根误差(rootmean squared error,RMSE＝0.1847mm)的重建精度。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。
Claims (9)
1.一种纳入重力的新型双平板X光机的校准方法，其特征在于，包括如下步骤：
S1：建立一套双平板X光机***，包括：
正位X光***、侧位X光***；
所述侧位X光***包括：第一平板X光机；
所述正位X光***包括：第二平板X光机；
所述第一平板X光机包括：第一X光源、第一平板探测器；
所述第二平板X光机包括：第二X光源、第二平板探测器；
所述第一平板X光机、第二平板X光机可固定地整体做上下移动；
S2：准备一校验盒；
所述校验盒包括依次相连的第一侧面、第二侧面、第三侧面、第四侧面；第一侧面和第三侧面对应第一平板X光机；第二侧面和第四侧面对应第二平板X光机；
所述第一侧面靠近第一平板探测器，第三侧面靠近第一X光源；将第一平板X光机拍摄的图片称为侧位图；
所述第二侧面靠近第二平板探测器，第四侧面靠近第二X光源；将第二平板X光机拍摄的图片称为正位图；
所述校验盒的每个侧面上均设置有若干个钢珠；且每个侧面的钢珠的分布方式一致；每个侧面的各个钢珠的三维坐标均已知；
S3：获取预存的校验盒的每个侧面各个钢珠的三维坐标；
S4:将校准盒摆在两台平板X光机光路重叠的区域，两台平板X光机同时曝光，获取正位图及侧位图；所述正位图和侧位图上均出现若干个黑点，每一个黑点均为一个对应钢珠的投影；
S5：进行侧位X光***和正位X光***校准前的对应关系准备；
S6：求出侧位X光***的内、外参数以及正位X光***的内、外参数；
S7：利用公式(1)求取正位X光***坐标系、侧位X光***坐标系的位姿关系；如下：
其中，记Rl,tl为侧位X光***的外参数，Rf,tf为正位X光***的外参数；
tl＝-RlCl；
tf＝-RfCf；
Pcf,Pcl是同一个点分别在正位X光***坐标系和侧位X光***坐标系下的坐标；
S8:将侧视图和正视图分别利用计算机图形学进行展示，即使用计算机图形学分别将侧位X光***、正位X光***的虚拟环境重现出来；
S9：确定重力方向，即确定重力向量在双平板X光***下的坐标表达。
2.如权利要求1所述的方法，其特征在于，所述步骤S5进一步包括：
S51：首先进行侧位X光***的校准准备工作，包括：
S511：使用“Hough圆检测”算法检测侧位图上黑点的圆心坐标，得到这些黑点的二维坐标，并将这些二维坐标集合在一起，记为二维点集xy，这里记录共有n个二维坐标；n为正整数；
S512：读取预存的钢珠坐标中涉及到侧位图的那一部分钢珠的坐标，将这些钢珠的坐标记为三维点集XYZ，这里记录共有n个三维坐标；
S513：建立从三维点集XYZ到二维点集xy的映射，得到二者的对应关系，共有n个对应关系；
S52：使用S51相同的方法，进行正位X光***的校准准备工作。
3.如权利要求2所述的方法，其特征在于，所述步骤S6进一步包括：
S61：首先求出侧位X光***的内、外参数，包括：
S611：建立空间三维点到图像二维点之间的变换矩阵：
从空间中三维点到图像上的二维点，借助齐次坐标，可用一次矩阵乘法完成：
其中：
Pw是步骤S51中三维点集XYZ的齐次形式，为4×1向量；p是像素坐标的齐次形式，为3×1向量；M是3×4矩阵，M矩阵蕴含了内参数和外参数；
S612：将第i个三维点集XYZ和二维点集xy的对应关系代入上式，并改写为齐次线性方程组的形式：
其中：
pi＝[ui vi]T是二维点集xy的第i个二维点；
Pi＝[Xi Yi Zi]T是三维点集XYZ中第i个三维点；
i＝1,…,n；
S613:将所有的n组对应关系联立，有如下过定线性齐次方程组：
记为
Am＝0
S614:对A奇异值分解：
A＝UDVT
U、D、V是A的三个分解矩阵；
过定方程的最优解是V矩阵的最后一列，即最优解
S615:构建非线性最优化问题：
其中：
pi＝[ui vi]T是二维点集xy的第i个二维点
Pi＝[Xi Yi Zi]T是三维点集XYZ中的第i个三维点坐标；
使用
作为此优化问题的初解，使用levernberg-marquardt最优化算法迭代求解，记优化终解为
S616:将向量
重排为矩阵形式，得到
S617:由
矩阵解析K,R,C；
其中：K是内矩阵，R是联系校准盒坐标系和侧位X光***坐标系的旋转矩阵，C是第一光源坐标；将K称为内参数，将R、C称为外参数；
将
矩阵左边3×3部分记为
则
对
使用QR分解，并将结果求逆可得K和R；
将
矩阵右边3×1部分记为
则
S62：使用S61相同的方法，求取正位X光***的内、外参数。
4.如权利要求1所述的方法，其特征在于，所述步骤S8进一步包括：
S81：将侧视图利用计算机图形学进行展示，即使用计算机图形学将侧位X光***的虚拟环境重现出来；进一步包括：
S811：将全局坐标系选在侧位X光***的原点；
S812：根据侧位X光***的内参数虚拟出一个虚拟成像平面，然后用纹理贴图的办法将实际获取的侧位图贴图到这个虚拟平面上；
S813：根据侧位X光***的外参数，将校验盒上涉及侧位成像的钢珠在空间中用黑色点显示出来；
S82：使用S81相同的方法，将正视图利用计算机图形学进行展示。
5.如权利要求1所述的方法，其特征在于，所述步骤S9进一步包括：
选取一重力方向的指征物，其同时被正位X光***、侧位X光***拍摄，分别获得所述指征物的正位图和侧位图；利用所述指征物在正位图和侧位图上的投影、以及之前求出的正位X光***、侧位X光***各自的内、外参数，推算重力方向。
6.如权利要求1所述的方法，其特征在于，所述步骤S9进一步包括：
选取一重力方向的指征物，并预先获得所述指征物的三维模型，使用正位X光***、侧位X光***对所述指征物进行拍摄，获得正位图、侧位图；
将预先获得的指征物的三维模型导入到正位X光***、侧位X光***的虚拟环境中；
将视点放在第二光源、第一光源处，同时手动调整所述指征物的三维模型在虚拟环境中的位置和朝向，当从视点处发现指征物的模型覆盖住了正位图、侧位图上指征物的轮廓，认为此时三维模型的位置和朝向信息还原出了指征物的实际位置和朝向信息；
根据指征物自身的形状特点、以及之前手动调整获得的指征物的实际空间位置和朝向信息，进一步推算重力方向信息。
7.如权利要求1所述的方法，其特征在于，还包括：
步骤S10：准备一验证盒，在S2的校验盒拍摄完毕后，保持第一平板X光机、第二平板X光机的位置不变，用验证盒替换校验盒，并拍摄正位图和侧位图；结合步骤S6获得的正、侧位X光***的内、外参数，使用“三角测量”技术，验证重建三维点的精度。
8.如权利要求7所述的方法，其特征在于，所述步骤S10进一步包括：
S101：使用单元块搭成一25cm*25cm*25cm的正方体盒子作为验证盒，将27颗钢珠嵌入正方体盒子内；
所述钢珠以近似3*3*3的晶格的形式分布于正方体内；钢珠之间在两个方向上间隔相同，在第三个方向间隔不同；
所述验证盒还内嵌入用钢丝制成的十字，圆圈，方形和三角标志；正、侧位图上的黑点即对应上述钢珠；
S102:预存所有27颗钢珠的三维坐标点，作为预存三维点集；
S103：使用“Hough”圆检测算法检测正、侧位图上的黑点；
将侧位片上的点集记为S100,将正位片上的点集记为S200；
S104:使用“三角测量”技术：于S100中的每一点，向第一光源引出一条射线，再取S200中的每一点，向第二光源引出一条射线，计算两射线的空间距离；
选择能产生最小射线空间距离的点对，认为这一对点是同一空间三维点在两侧图像上的投影；
同时计算两射线最靠近时，最小二乘意义上的交点，作为重建点三维坐标；
循环使用此方法，可建立S100与S200的对应关系，对每一点对，也同时重建了其对应的三维点坐标；利用以上方法计算出来的27颗钢珠三维点坐标的集合，即为重建的三维点集；
S105:步骤S104的重建的三维点集是位于双平板X光***坐标系下的，步骤S102的预存三维点集是位于验证盒坐标系下的，故将步骤S102的预存点集和步骤S104的重建的三维点集配准到同一坐标系下，方法如下：
结合步骤S6获得的正、侧位X光***的内、外参数，对重建三维点集使用PCA确定其中心，从而确定两个坐标系的偏移向量；
PCA给出的前三个主成分分量其实就是旋转过的三坐标轴；
利用钢珠排列在两个方向上间隔相同在第三个方向间隔不同这一空间异质性，配合可视化来看是否需要坐标轴的翻转来更好的拟合两组点；
当从可视化中看出重建的三维点集和预存的三维点集基本重合后，再用暴力求解法确定点的对应关系；
最终可计算出RMSE来表征校准精度。
9.如权利要求1所述的方法，其特征在于，
第一平板X光机、第二平板X光机的光轴存在一夹角，所述夹角为30度-150度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010775024.XA CN111973204B (zh) | 2020-08-04 | 2020-08-04 | 一种纳入重力的新型双平板x光机的校准方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010775024.XA CN111973204B (zh) | 2020-08-04 | 2020-08-04 | 一种纳入重力的新型双平板x光机的校准方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111973204A CN111973204A (zh) | 2020-11-24 |
| CN111973204B true CN111973204B (zh) | 2022-09-06 |
Family
ID=73444971
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010775024.XA Active CN111973204B (zh) | 2020-08-04 | 2020-08-04 | 一种纳入重力的新型双平板x光机的校准方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111973204B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112419381B (zh) * | 2020-12-15 | 2023-03-03 | 山东威高医疗科技有限公司 | 一种x光图像中标记点序列的自动识别方法 |
| CN114511641A (zh) * | 2021-01-30 | 2022-05-17 | 威海威高骨科手术机器人有限公司 | 带有位置信息的x光数据二维标定处理方法 |
| CN112927296A (zh) * | 2021-02-03 | 2021-06-08 | 上海橙捷健康科技有限公司 | 一种用于空间相对位置标定与校准的方法及*** |
| CN113643428A (zh) * | 2021-08-17 | 2021-11-12 | 北京唯迈医疗设备有限公司 | 一种适用于多自由度锥形束ct的全参数几何校准方法 |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101149836B (zh) * | 2007-11-05 | 2010-05-19 | 中山大学 | 一种三维重构的双摄像机标定方法 |
| KR101107513B1 (ko) * | 2009-10-14 | 2012-01-31 | (주) 케이앤아이테크놀로지 | 이중 방사선 촬영 장비의 캘리브레이션 방법 및 이를 이용한 수술 후 관절 의 삼차원 위치 정보의 획득 방법 |
| CN102743184B (zh) * | 2012-05-14 | 2013-10-16 | 清华大学 | 一种x射线锥束计算机层析成像***的几何参数标定方法 |
| CN104581136B (zh) * | 2013-10-14 | 2017-05-31 | 钰立微电子股份有限公司 | 图像校准***和立体摄像机的校准方法 |
| CN105606127A (zh) * | 2016-01-11 | 2016-05-25 | 北京邮电大学 | 一种双目立体相机与惯性测量单元相对姿态标定方法 |
| CN107133989B (zh) * | 2017-06-12 | 2020-11-06 | 中国科学院长春光学精密机械与物理研究所 | 一种三维扫描***参数标定方法 |
| CN108154538A (zh) * | 2018-02-06 | 2018-06-12 | 华中科技大学 | 一种双摄像机模组校正和标定方法及装置 |
-
2020
- 2020-08-04 CN CN202010775024.XA patent/CN111973204B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN111973204A (zh) | 2020-11-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111973204B (zh) | 一种纳入重力的新型双平板x光机的校准方法 | |
| US20240123260A1 (en) | Method of calibration of a stereoscopic camera system for use with a radio therapy treatment apparatus | |
| CN104783824B (zh) | X射线成像***的校正方法 | |
| EP1278458B1 (en) | Fluoroscopic tracking and visualization system | |
| DE102006005036B4 (de) | Vorrichtung und Verfahren zur verbesserten Formcharakterisierung | |
| EP3453330A1 (en) | Virtual positioning image for use in imaging | |
| EP3682424B1 (en) | Relative pose determination of non-overlapping cameras | |
| US20050151963A1 (en) | Transprojection of geometry data | |
| Zhang et al. | A universal and flexible theodolite-camera system for making accurate measurements over large volumes | |
| CN106456100A (zh) | 用于配置x射线成像***的方法和*** | |
| CA2926529A1 (en) | Tracking system for imaging machines and related apparatus | |
| CN112168357B (zh) | C臂机空间定位模型构建***及方法 | |
| US10687900B2 (en) | Method of image support for a person carrying out a minimally invasive procedure with an instrument in a procedure site of a patient, X-ray apparatus, computer program and electronically readable data carrier | |
| Selby et al. | Patient positioning with X-ray detector self-calibration for image guided therapy | |
| CN110910506B (zh) | 基于法线检测的三维重建方法、装置、检测装置及*** | |
| CN109953767A (zh) | 用于校准医学成像设备、用于执行配准的方法以及*** | |
| CN113298886A (zh) | 一种投影仪的标定方法 | |
| KR101524937B1 (ko) | 방사선 캘리브레이션 방법 및 그에 의한 캘리브레이션 장치 | |
| US7241045B2 (en) | Stereoradiography device and method for the use thereof | |
| KR102479266B1 (ko) | 치료 시스템, 캘리브레이션 방법, 및 프로그램 | |
| Cheon et al. | High-precision quality assurance of robotic couches with six degrees of freedom | |
| CN113963057B (zh) | 成像几何关系标定方法、装置、电子设备以及存储介质 | |
| CN113963058B (zh) | 预定轨道ct在线标定方法、装置、电子设备及存储介质 | |
| Cheryauka et al. | 3-D geometry calibration and markerless electromagnetic tracking with a mobile C-arm | |
| Bujakiewicz et al. | Modelling and visualization of three dimensional objects using close range imagery |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20211223 Address after: 201210 room 405-2, 4th floor, building 9, No. 1206, Zhangjiang Road, China (Shanghai) pilot Free Trade Zone, Pudong New Area, Shanghai Applicant after: SHANGHAI TAOYING MEDICAL TECHNOLOGY CO.,LTD. Address before: 200240 No. 800, Dongchuan Road, Shanghai, Minhang District Applicant before: SHANGHAI JIAO TONG University |
|
| TA01 | Transfer of patent application right | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |