CN104764756A - 锥束ct成像***的标定方法 - Google Patents

Abstract

本发明公开了一种锥束CT成像***的标定方法，包括：建造模体，模体大体上是塑料圆柱体，模体的上下底面各镶嵌了两组钢珠，每组包括N个钢珠；获取2N个钢珠中心的投影坐标，并分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂；得到i系原点O和交点D并求出角度η、A₁、B₁、A₂、B₂四点坐标；构造第一目标函数f(θ*,w_y*,w_z*)以计算角度θ和在i系下，计算W坐标系的原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)；构造第二目标函数并计算角度构造第三目标函数H(t)并求出角度t；计算在W系下的光源P_s坐标和原点O坐标()。本发明能对任意复杂轨道的CBCT***进行标定，且无需预先知道***的任何几何信息，利用一次投影即可得到该投影位置的全部几何参数，普适性和针对性强。

Description

锥束CT成像***的标定方法

技术领域

本发明涉及X光计算机断层成像技术领域，尤其涉及一种锥束CT成像***的标定方法。

背景技术

锥束CT(Cone Beam Computed Tomography，CBCT)技术是应用锥束X光源围绕三维物体作各个角度的旋转，获取物体在多个角度下的投影，然后运用锥束投影图的重建算法重建出物体三维图像的技术。其原理主要是利用物体不同成分对X光的衰减系数不同的特点，在物体后方探测器上得到能反映物体组分特征的投影图像。

CBCT***的几何参数是指描述***的X光源、物体和探测器之间的角度和位置的几何参量。这些几何参量会影响物体的投影图像，并会作为重要的信息在重建算法中使用。因此CBCT***的几何参数标定具有提高重建图像质量，防止出现图像伪影的重要意义。

在实际应用的过程中，CBCT***中X光源的扫描轨迹不一定是完整的圆周，比如相关技术中C-arm CT***的光源被固定在C型旋转臂上，医生可以在为患者手术的同时，及时观察成像情况。这种C-arm CT***的扫描轨迹不是一个完整的圆周，另外由于重力的作用，轨迹的中心会出现偏离圆心的情况，造成X光源到探测器的距离发生改变，这种沿不规则轨道扫描的CBCT***同样需要精确标定它的几何参数。

发明内容

本发明旨在至少在一定程度上解决相关技术中的技术问题之一。为此，本发明的一个目的在于提出一种锥束CT成像***的标定方法，该锥束CT成像***的标定方法可以对任意复杂轨道的CBCT***进行标定，具有较强的普适性。

为了实现上述目的，本发明第一方面实施例的锥束CT成像***的标定方法，包括以下步骤：建造模体，其中，所述模体为标定参照物，所述模体大体上是塑料圆柱体，所述模体的上下底面各镶嵌了两组钢珠，每组包括N(N不小于5)个钢珠，所述模体底面圆的直径记为2r，所述两组钢珠中心所决定的两个平面的间距记为2l；以所述钢珠投影在系的所有像素的值为权重对每个钢珠的投影求加权平均坐标值，以获取所述2N个钢珠中心的投影坐标，并分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂；分别将所述每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序，并得到系下的原点O和交点D以及计算角度η和A₁、B₁、A₂、B₂四点坐标，其中，所述2N个钢珠中心投影坐标均在所述系下表达，所述A₁、B₁分别为P₁、Q₁在OD上的投影，所述A₂、B₂分别为P₂、Q₂在所述OD上的投影，称矩形P₁P₂Q₁Q₂(长为2l，宽为2r)为探测器绕所述系下的旋转角度θ后，所述模体被系下的平面所截得的截面，所述角度η为所述探测器在绕所述旋转所述角度θ和绕所述系下的旋转角度后，绕所述系下的旋转的角度；在所述平面内，根据所述角度η、所述原点O、所述交点D、Q₁坐标(s₁,t₁)、P₂坐标(s₂,t₂)、A₁坐标A₂坐标B₁坐标B₂坐标构造第一目标函数f(θ*,w_y*,w_z*)，并以此计算所述角度θ、在所述平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)，其中，所述原点O_w同时也是所述模体的中心；根据所述模体两个底面圆上任意一点的投影应分别位于所述轨迹T₁和所述轨迹T₂上的几何关系来构造第二目标函数并计算所述角度根据所述原点O、所述角度θ、所述角度η、所述角度和在所述平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)以及钢珠投影中心坐标间的位置关系构造第三目标函数H(t)并求出角度t；以及根据所述原点O、所述角度θ、所述角度η、所述角度所述角度t和在所述平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)计算在系下的光源P_s坐标和原点O坐标

根据本发明实施例的锥束CT成像***的标定方法，具有以下有益效果：(1)可以对任意复杂轨道的CBCT***进行标定，具有较强的普适性。(2)利用一次投影即可得到该投影位置的全部几何参数，从而对于CT***扫描角度不全的情况非常有针对性。(3)不需要预先知道任何***的几何信息，可以用于自身机械精度很差的几何***。

进一步地，在本发明的一个实施例中，所述系为相对周围环境静止的世界坐标系，所述系以所述模体的中心为原点O_w，所述系以所述模体的轴为在开始投影前，所述系以垂直于所述的射线为所述系下的根据所述所述和右手定则确定；所述系为描述所述探测器理想状态的虚拟探测器坐标系，所述系以所述光源P_s和所述原点O_w的连线与所述探测器平面的交点为所述原点O，所述和所述同向，定义P_sO与所述确定的平面为所述平面，所述过所述原点O且垂直于所述根据所述所述和右手定则确定；以及所述系为基于所述探测器本身的坐标系，所述系以所述原点O为原点，所述和所述平行于所述探测器的行列线，所述垂直于所述探测器平面，其中，定义所述所述所述分别与所述所述所述对应重合时的状态为所述探测器的初始状态。

进一步地，在本发明的一个实施例中，从所述探测器的初始状态到所述探测器的任一状态包括以下步骤：所述探测器从所述探测器的初始状态绕所述旋转，与所述形成所述角度θ，此处，所述角度θ值按照以所述为正方向的右手定则方式定义；所述探测器从所述探测器的初始状态绕所述旋转，与所述形成所述角度此处，所述角度值按照以所述为正方向的右手定则方式定义；以及所述探测器从所述探测器的初始状态绕所述旋转，与所述形成所述角度η，此处，所述角度η值按照以所述为正方向的右手定则方式定义。

进一步地，在本发明的一个实施例中，根据以下方程分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂：

$a_1{(u-u_{0_1})}^2+b_1{(v-v_{0_1})}^2+2c_1(u-u_{0_1})(v-v_{0_1})-1=0$

$a_2{(u-u_{0_2})}^2+b_2{(v-v_{0_2})}^2+2c_2(u-u_{0_2})(v-v_{0_2})-1=0$

其中，a₁、b₁、c₁、和a₂、b₂、c₂、由每组的所述N个钢珠中心投影在所述系的坐标拟合得到。

进一步地，在本发明的一个实施例中，所述分别将所述每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序，并得到系下的原点O和交点D具体为：分别将所述每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序；连接E₁F₂和E₂F₁得到所述原点O，其中，所述E₁和所述E₂分别为所述轨迹T₁上纵坐标最大的钢珠中心投影和纵坐标最小的钢珠中心投影，所述F₁和所述F₂分别为所述轨迹T₂上纵坐标最大的钢珠中心投影和纵坐标最小的钢珠中心投影；以及连接E₁F₁和E₂F₂，所述E₁F₁和所述E₂F₂在空间中延长线的交点作为所述交点D，其中，所述OD与所述的夹角为所述角度η。

进一步地，在本发明的一个实施例中，根据以下步骤计算所述角度θ、在所述平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)：在所述平面内，对于待定的原点O_w*坐标(w_y*,w_z*)，根据所述矩形P₁P₂Q₁Q₂的长和宽得到相应的待定矩形顶点P₁*、P₂*、Q₁*、Q₂*，连接所述原点O和所述原点O_w*得到直线OO_w*；待定角度θ*，过所述交点D作所述的平行线DP_s*和所述直线OO_w*相交于点P_s*，连接A₁P₁*、B₁Q₁*、A₂P₂*、B₂Q₂*分别与所述直线OO_w*相交于构造第一目标函数f(θ*,w_y*,w_z*)用以表示所述P_s*、的离散程度，其中，所述第一目标函数f(θ*,w_y*,w_z*)越大，所述P_s*、的离散程度越大，当所述第一目标函数f(θ*,w_y*,w_z*)＝0时表示所述P_s*、五点重合；计算 $({\mathrm{\θ}}_0,w_{y_0},w_{z_0})=\arg \min \left\{f({\mathrm{\θ}}^{*},{w_y}^{*},{w_z}^{*})\right\}$ 的优化结果，以获得所述角度θ和在所述平面内的所述原点O_w坐标(w_y,w_z)，即：θ＝θ₀；取θ₀,对应求出的所述P_s*、 五点的平均坐标为所述光源P_s坐标(s_y,s_z)。

进一步地，在本发明的一个实施例中，或 $f(\mathrm{\θ}*,w_y*,w_z*)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_s*P_{s_i}\left|\right|}^{\mathrm{\λ}}+\underset{i,j=1,i\≠j}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_{s_i}{P}_{s_j}\left|\right|}^{\mathrm{\λ}}$ 其中，λ一般为1或2。

进一步地，在本发明的一个实施例中，根据以下步骤计算所述角度待定角度根据已求出的所述角度η、所述角度θ、在所述平面内的所述原点O_w坐标(w_y,w_z)、所述光源P_s坐标(s_y,s_z)和待定角度计算所述模体两个底面圆在所述系的轨迹T₁*、T₂*；在所述轨迹T₁*、T₂*上，分别取点(u₁,v₁)和(u₂,v₂)，构造函数以衡量所述点(u₁,v₁)和(u₂,v₂)分别与所述轨迹T₁和所述轨迹T₂的偏离程度，其中，所述越大，所述点(u₁,v₁)和(u₂,v₂)分别与所述轨迹T₁和轨迹T₂的偏离程度越大；构造第二目标函数以计算所述角度其中，M表示共取了M组在所述轨迹T₁*、T₂*上的所述点(u₁,v₁)和(u₂,v₂)，

进一步地，在本发明的一个实施例中，

其中，a₁、b₁、c₁、和a₂、b₂、c₂、由每组的所述N个钢珠中心投影在所述系的坐标拟合得到。

进一步地，在本发明的一个实施例中，根据以下公式计算所述角度t：

t＝argmin{H(t*)}

其中，角度参数t*描述所述模体底面圆上的钢珠位于该底面圆上的位置，函数h₁(t*)、h₂(t*)分别为在已知所述角度η、所述角度θ、所述原点O_w坐标(w_y,w_z)，所述光源P_s坐标(s_y,s_z)和所述角度的前提下，在模体两个底面上处于t*所对应位置的钢珠的投影和在所述系下最近的钢珠投影中心的距离。

进一步地，在本发明的一个实施例中，

附图说明

图1是根据本发明一个实施例的锥束CT成像***的标定方法的流程图；

图2是根据本发明一个实施例的锥束CT成像***的标定方法中模体的结构示意图；

图3是根据本发明一个实施例的锥束CT成像***的标定方法的坐标系定义示意图；

图4是根据本发明一个实施例的锥束CT成像***的标定方法中探测器从探测器的初始状态绕旋转角度θ的旋转示意图；

图5是根据本发明一个实施例的锥束CT成像***的中探测器继续绕旋转角度的旋转示意图；

图6是根据本发明一个实施例的锥束CT成像***的标定方法中探测器继续绕旋转角度η的旋转示意图；

图7是根据本发明一个实施例的锥束CT成像***的标定方法中探测器继续绕旋转角度η后系和系的坐标系示意图；

图8是根据本发明一个实施例的锥束CT成像***的标定方法中的平面的示意图；

图9是根据本发明一个实施例的锥束CT成像***的标定方法的原点O和交点D的定位示意图；以及

图10是根据本发明一个具体实施例的锥束CT成像***的标定方法的流程。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

下面参考附图描述本发明实施例的锥束CT成像***的标定方法。

图1是根据本发明一个实施例的锥束CT成像***的标定方法的流程图。如图1所示，本发明实施例的锥束CT成像***的标定方法，包括以下步骤：

S1，建造模体，其中，模体为标定参照物，模体大体上是塑料圆柱体，模体的上下底面各镶嵌了两组钢珠，每组包括N(N不小于5)个钢珠，模体底面圆的直径记为2r，两组钢珠中心所决定的两个平面的间距记为2l。

具体地，在本发明的一个实施例中，如图2所示，两组钢珠中心分别设置在模体的两底面圆上，且每组的N个(例如6个)钢珠等角均匀放置，两组钢珠中心所决定两底面圆的平面间距为2l，每组的N个钢珠中心决定的底面圆的直径为2r，N、l、r和钢珠自身的半径需满足能够在投影图上识别钢珠中心位置的前提下，使得投影面积尽可能大。

进一步地，对于任意待标定的CBCT***，首先定义如下3个坐标系系、系和系。如图3所示，系为相对周围环境静止的世界坐标系，系以模体的中心为原点O_w，系以模体的轴为(y轴，的正反向任意)，在开始投影前，系以垂直于的某条射线为(x轴)，系下的(z轴)根据和右手定则确定。系为描述探测器理想状态的虚拟探测器坐标系，系以光源P_s和原点O_w的连线与探测器平面的交点为原点O，(y轴)和同向，定义P_sO与确定的平面为平面，(z轴)过原点O且垂直于(x轴)根据和右手定则确定，其中，的正方向为光源P_s在上的投影位于正半轴时的方向。系为基于探测器本身的坐标系，系以原点O为原点，(x轴)和(y轴)平行于探测器的行列线，(z轴)垂直于探测器平面。其中，定义分别与对应重合时的状态为探测器的初始状态。

进一步地，在本发明的一个实施例中，从探测器的初始状态到探测器的任一状态包括以下步骤：

S8，探测器从探测器的初始状态绕旋转，与形成角度θ，此处，角度θ值按照以为正方向的右手定则方式定义。

步骤S8中探测器的旋转示意图如图4所示。

S9，探测器从探测器的初始状态绕旋转，与形成角度此处，角度值按照以为正方向的右手定则方式定义。

步骤S9中探测器的旋转示意图如图5所示。

S10，探测器从探测器的初始状态绕旋转，与形成角度η，此处，角度η值按照以为正方向的右手定则方式定义。

步骤S10中探测器的旋转示意图如图6所示，此时的系和系如图7所示。

需要说明的是，步骤S8至步骤S10中，探测器的每一次旋转都会使得中两个轴的方向发生变化，为使图4至图6更具可理解性，图4至图6中是指探测器该次旋转之前的也即旋转角度θ、角度角度η所分别围绕的

S2，以钢珠投影在系的所有像素的值为权重对每个钢珠的投影求加权平均坐标值，以获取2N个钢珠中心的投影坐标，并分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂。

由于每组N个钢珠在空间中的位置处于底面圆上，所以每组的N个钢珠中心投影坐标应该处于底面圆的投影上，即一个椭圆上。在本发明的一个实施例中，可以根据以下方程分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂：

$a_1{(u-u_{0_1})}^2+b_1{(v-v_{0_1})}^2+2c_1(u-u_{0_1})(v-v_{0_1})-1=0$

$a_2{(u-u_{0_2})}^2+b_2{(v-v_{0_2})}^2+2c_2(u-u_{0_2})(v-v_{0_2})-1=0$

其中，a₁、b₁、c₁、和a₂、b₂、c₂、由每组的N个钢珠中心投影在系的坐标拟合得到。

S3，分别将每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序，并得到系下的原点O和交点D以及计算角度η和A₁、B₁、A₂、B₂四点坐标，其中，2N个钢珠中心投影坐标均在系下表达，A₁、B₁分别为P₁、Q₁在OD上的投影，A₂、B₂分别为P₂、Q₂在OD上的投影，称矩形P₁P₂Q₁Q₂(长为2l，宽为2r)为探测器绕系下的旋转角度θ后，模体被系下的平面所截得的截面，角度η为探测器在绕旋转角度θ和绕系下的旋转角度后，绕系下的旋转的角度。

需要说明的是，如图8所示，探测器绕系下的旋转角度θ后，和OD所在直线重合且和的夹角为角度θ，探测器绕旋转角度θ和绕旋转角度并绕系下的轴旋转角度η后，OD和的夹角为角度η。图8中，||P₁P₂||＝2l，||P₁Q₁||＝2r，直线P_sD平行于并与探测器平面相交于D，A₁B₁、A₂B₂分别由轨迹T₁和轨迹T₂与OD相交得到。

进一步地，在本发明的一个实施例中，分别将每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序，并得到系下的原点O和交点D具体可以为：

S31，分别将每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序。

S32，连接E₁F₂和E₂F₁得到原点O，其中，E₁和E₂分别为轨迹T₁上纵坐标最大的钢珠中心投影和纵坐标最小的钢珠中心投影，F₁和F₂分别为轨迹T₂上纵坐标最大的钢珠中心投影和纵坐标最小的钢珠中心投影。

E₁、E₂、F₁和F₂如图3所示。

S33，连接E₁F₁和E₂F₂，E₁F₁和E₂F₂在空间中延长线的交点作为交点D，其中，OD与的夹角为角度η。

步骤S32和步骤S33中定位原点O和交点D的示意图如图9所示。

S4，在平面内，根据角度η、原点O、交点D、Q₁坐标(s₁,t₁)、P₂坐标(s₂,t₂)、A₁坐标A₂坐标B₁坐标B₂坐标构造第一目标函数f(θ*,w_y*,w_z*)，并以此计算角度θ、在平面内的原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)，其中，原点O_w同时也是模体的中心。

具体地，在本发明的一个实施例中，可以根据以下步骤计算角度θ、在平面内的原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)：

S41，在平面内，对于待定的原点O_w*坐标(w_y*,w_z*)，根据矩形P₁P₂Q₁Q₂的长和宽得到相应的待定矩形顶点P₁*、P₂*、Q₁*、Q₂*，连接原点O和原点O_w*得到直线OO_w*。

S42，待定角度θ*，过交点D作的平行线DP_s*和直线OO_w*相交于点P_s*，连接A₁P₁*、B₁Q₁*、A₂P₂*、B₂Q₂*分别与直线OO_w*相交于

S43，构造第一目标函数f(θ*,w_y*,w_z*)用以表示P_s*、的离散程度，其中，第一目标函数f(θ*,w_y*,w_z*)越大，P_s*、的离散程度越大，当第一目标函数f(θ*,w_y*,w_z*)＝0时表示P_s*、五点重合。

S44，计算 $({\mathrm{\θ}}_0,w_{y_0},w_{z_0})=\arg \min \left\{f({\mathrm{\θ}}^{*},{w_y}^{*},{w_z}^{*})\right\}$ 的优化结果，以获得角度θ和在平面内的原点O_w坐标(w_y,w_z)，即：θ＝θ₀。

S45，取θ₀,对应求出的P_s*、五点的平均坐标为光源P_s坐标(s_y,s_z)。

具体地，在本发明的一个实施例中， $(y_{A_1},z_{A_1})=\left(\right|\left|O{A}_1\right||\cos (\mathrm{\θ}+\mathrm{\π}),|\left|O{A}_1\right|\left|\sin (\mathrm{\θ}+\mathrm{\π})\right)$ $(y_{B_1},z_{B_1})=\left(\right|\left|O{B}_1\right||\cos (\mathrm{\θ}+\mathrm{\π}),|\left|O{B}_1\right|\left|\sin (\mathrm{\θ}+\mathrm{\π})\right),(y_{A_2},z_{A_2})=\left(\right|\left|O{A}_2\right||\cos \mathrm{\θ},|\left|O{A}_2\right|\left|\sin \mathrm{\θ}\right)$ $(y_{B_2},z_{B_2})=\left(\right|\left|O{B}_2\right||\cos \left(\mathrm{\θ}\right),|\left|O{B}_2\right|\left|\sin \left(\mathrm{\θ}\right)\right),s_y=\frac{w_y{s}_z}{w_z},$ s_z＝||OD||sinθ，(s₁,t₁)＝(w_y-l,w_z-r)，(s₂,t₂)＝(w_y+l,w_z+r)。进一步地，根据P_sQ₁B₁，P_sP₂A₂三点共线，可知： $(z_{B_1}-s_z)(y_{B_1}-s_1)=(z_{B_1}-t_1)(y_{B_1}-s_y),(z_{A_2}-s_z)(y_{A_2}-s_2)=(z_{A_2}-t_2)(y_{A_2}-s_y).$ 令 $k=\frac{y_{B_1}}{y_{A_2}}=\frac{z_{B_1}}{z_{A_2}}.$

进一步地，在本发明的一个实施例中，当时，O为B₁A₂中点，B₁A₂平行Q₁P₂，其中，2r为底面圆的直径，2l为两组钢珠中心所决定的两个平面的间距。

进一步地，在本发明的另一个实施例中，当时， $\mathrm{\ξ}=\frac{w_y}{w_z}=\frac{C-\mathrm{kD}}{r\left|\right|\mathrm{OD}\left|\right|\sin \mathrm{\θ}(1+k)},w_z=\frac{C-\mathrm{\ξ}\left|\right|\mathrm{OD}\left|\right|\sin \mathrm{\θ}(z_{B_1}+r)}{y_{B_1}-\mathrm{\ξ}{z}_{B_1}},$ w_y＝ξw_z， $\mathrm{\θ}=\underset{\mathrm{\θ}}{\arg \min }\left\{f(\mathrm{\θ}*,w_y*,w_z*)\right\},$ 其中， $C=r{y}_{B_1}+\left|\right|\mathrm{OD}\left|\right|y_{B_1}\sin \mathrm{\θ}-z_{B_1}l+\left|\right|\mathrm{OD}\left|\right|l\sin \mathrm{\θ},$ $D=r{y}_{A_2}+\left|\right|\mathrm{OD}\left|\right|y_{A_2}\sin \mathrm{\θ}-z_{A_2}l+\left|\right|\mathrm{OD}\left|\right|l\sin \mathrm{\θ},f(\mathrm{\θ}*,w_y*,w_z*)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_s*P_{s_i}\left|\right|}^{\mathrm{\λ}},$ 或 $f(\mathrm{\θ}*,w_y*,w_z*)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_s*P_{s_i}\left|\right|}^{\mathrm{\λ}}+\underset{i,j=1,i\≠j}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_{s_i}*P_{s_j}\left|\right|}^{\mathrm{\λ}}$ 等，其中，λ一般为1或2。

S5，根据模体两个底面圆上任意一点的投影应分别位于轨迹T₁和轨迹T₂上的几何关系来构造第二目标函数并计算角度

进一步地，在本发明的一个实施例中，可以根据以下步骤计算角度

S51，待定角度根据已求出的角度η、角度θ、在平面内的原点O_w坐标(w_y,w_z)、光源P_s坐标(s_y,s_z)和待定角度计算模体两个底面圆在系的轨迹T₁*、T₂*。

S52，在轨迹T₁*、T₂*上，分别取点(u₁,v₁)和(u₂,v₂)，构造函数以衡量点(u₁,v₁)和(u₂,v₂)分别与轨迹T₁和轨迹T₂的偏离程度，其中，越大，点(u₁,v₁)和(u₂,v₂)分别与轨迹T₁和轨迹T₂的偏离程度越大。

具体地，在本发明的一个实施例中，

其中，a₁、b₁、c₁、和a₂、b₂、c₂、由每组的N个钢珠中心投影在系的坐标拟合得到。

S53，构造第二目标函数以计算角度其中，M表示共取了M组在轨迹T₁*、T₂*上的点(u₁,v₁)和(u₂,v₂)，

其中，M为大于或等于1的整数，j为小于M的整数，第二目标函数用于衡量角度下钢珠中心投影坐标与轨迹T₁和轨迹T₂的偏差。

需要说明的是，平面绕旋转-η时，平面上所有点的位置也都绕旋转了-η。此时，在系下，新的系中三个基向量分别可以表示为：

进一步地，在本发明的一个实施例中，对于两底面圆上的分别任意一点K_1,2，在系下，其中，α_1,2∈[0,360°)，定义与平面的交点为R_1,2，则在系下，且解得：根据λ求出在系下的坐标，再用即可得到R_1,2在新的系下的坐标将坐标绕原点O旋转角度η后再分别代入拟合轨迹T₁和轨迹T₂的方程，因此，根据轨迹T₁和轨迹T₂上的任意多个点即可获得轨迹T₁和轨迹T₂的方程中各参数。

S6，根据原点O、角度θ、角度η、角度和在平面内的原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)以及钢珠投影中心坐标间的位置关系构造第三目标函数H(t)并求出角度t。

其中，角度t为系下的正半轴在系下的平面上的投影和系下的正方向的夹角，且角度t的正方向和系下的满足右手定则，其中，模体的第一张投影中和重合。

由于模体的第一张投影对应***的和重合，所以模体的第一张投影对应***的t＝0。进一步地，在本发明的一个实施例中，可以根据以下公式计算角度t：

$t=\underset{t}{\arg \min }\left\{H\left(t\right)\right\}$

其中，角度参数t*描述模体底面圆上的钢珠位于该底面圆上的位置，函数h₁(t*)、h₂(t*)分别为在已知角度η、角度θ、原点O_w坐标(w_y,w_z)，光源P_s坐标(s_y,s_z)和角度的前提下，在模体两个底面上处于t*所对应位置的钢珠的投影和在系下最近的钢珠投影中心的距离。具体地，在本发明的一个实施例中， $h_1\left(t\right)=\sqrt{{(x_{R_1}-x_{J_1})}^2+{(y_{R_1}-y_{J_1})}^2}$ 为绕原点O旋转角度η后，与相距最近的钢珠中心的距离， $h_2\left(t\right)=\sqrt{{(x_{R_2}-x_{J_2})}^2+{(y_{R_2}-y_{J_2})}^2}$ 为绕原点O旋转角度η后，与相距最近的钢珠中心的距离，为P_sK_1,2与平面的交点，K_1,2为两底面圆上的分别任意一点，在系下，i为大于等于0且小于等于N-1的整数。H(t)以为周期，确定H(t)的真值需要结合前一次投影标定的角度t。

S7，根据原点O、角度θ、角度η、角度角度t和在平面内的原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)计算在系下的光源P_s坐标和原点O坐标

进一步地，在本发明的一个实施例中，系下的光源P_s坐标和原点O的坐标可以根据以下公式进行计算：

以下根据具体实施例对本发明实施例的锥束CT成像***的标定方法进行说明。

如图10所示，在本发明的一个具体实施例中，锥束CT成像***的标定方法包括以下步骤：

S101，建造模体。

其中，测试模体取N＝6，l＝40mm，r＝65mm，钢珠自身的直径为2.5mm，钢珠总数目为12。待标定***的探测器为1920×1536像素，探测器像素单元边长0.127mm。对于任意待标定的CBCT***，定义坐标系系、系、系，探测器的旋转角θ、η，以及O到垂线绕的旋转角度t。

S102，以钢珠投影所有像素的值为权重对每个钢珠的投影求加权平均坐标值，得到12个钢珠的球心坐标，拟合出12个钢珠的球心坐标所在椭圆的轨迹T₁，T₂。

S103，将每组的6个钢珠中心定好顺序，并得到原点O，此时，所有12个钢珠中心坐标都可以在系下表达。

S104，计算出角度η，并由轨迹T₁、T₂和OD相交得到A₁B₁、A₂B₂。

S105，在平面内，求出角度θ、角度θ对应的原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)。

角度θ误差在0.01°左右，但是关于原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)的精度并不理想(尤其在θ→0,OD→∞时)。

S106，对角度θ、原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)作第一次标定修正。

如果角度θ>2°(2°由实验条件确定)，则原点的坐标就取已求得的(w_y,w_z)。反之作如下处理：把角度θ近似为0。假设P₁Q₁中点为M₁，P₂Q₂中点为M₂，则根据“两直线交点的投影等于其投影的交点”，可以在椭圆T₁、T₂中用连接相对钢珠中心并取连线交点的方法求出M₁、M₂的投影O₁、O₂(O₁、O₂在直线OD上)。

令||OP_s||/||OO_w||＝p'，即(s_y,s_z)＝p'(w_y,w_z)，则||O₁O₂||/||M₁M₂||＝p'/(p'-1)，得：

$p^{\text{'}}=\frac{1}{\left|\right|O_1{O}_2\left|\right|/2l-1}+1,s_z=\frac{r}{1-\frac{1}{p^{\text{'}}}-\frac{2l}{\left|\right|O{A}_1\left|\right|+\left|\right|O{A}_2\left|\right|}},$

$s_y=\frac{\left(\right|\left|O{A}_1\right||-|\left|O{A}_2\right|\left|\right)(s_z-\frac{s_z}{p^{\text{'}}}-r)+\left(\right|\left|O{B}_2\right||-|\left|O{B}_1\right|\left|\right)(s_z-\frac{s_z}{p^{\text{'}}}+r)}{4r}(w_y,w_z)=\frac{1}{p^{\text{'}}}(s_y,s_z).$

S107，对原点O_w坐标(w_y,w_z)和光源P_s坐标(s_y,s_z)进行进一步的标定修正。

步骤S107也可以重复多次，具体方法是：

在(w_y±0.2,w_z±0.2)(0.2mm为该实验条件下第一次标定修正后的最大误差)的范围之内取求出此时的P₁P₂Q₁Q₂的坐标(w_y'±l,w_z'±r)，然后连接A₁P₁、B₁Q₁、A₂P₂、B₂Q₂与相交于目标函数为 $g({w_y}^{\text{'}},{w_z}^{\text{'}})=\mathrm{var}\left(\right|\left|O{P}_{s{A}_1}\right||,|\left|O{P}_{s{B}_1}\right||,|\left|O{P}_{s{A}_2}\right||,|\left|O{P}_{s{B}_2}\right|\left|\right)$ 表示求的方差)，取 $w_y,w_z=\underset{w_y,w_z}{\arg \min }\left\{g(w_y,w_z)\right\},$ w_y,w_z所对应的 $(s_y,s_z)=\frac{(s_{y{A}_1},s_{z{A}_1})+(s_{y{B}_1},s_{z{B}_1})+(s_{y{A}_2},s_{z{A}_2})+(s_{y{B}_2},s_{z{B}_2})}{4}$ 即为要求的(s_y,s_z)。

S108，得到新的系，并定义子函数定义目标函数来衡量角度下球心投影和T₁,T₂轨迹的偏差，计算角度

取：这时α₁，α₂＝10j，j＝0,1,2…35在椭圆T₁,T₂上。

S109，定义实验条件下的目标函数H(t)计算角度t。

$H\left(t\right)=\underset{i=0}{\overset{5}{\mathrm{\Σ}}}{h}_1(t+60i)+h_2(t+60i),$ 取： $t=\underset{t}{\arg \min }\left\{H\left(t\right)\right\}.$

S110，求出在系下光源P_s坐标和原点O坐标

S111，标定结束。

本发明的有益效果如下：(1)可以对任意复杂轨道的CBCT***进行标定，具有较强的普适性。(2)利用一次投影即可得到该投影位置的全部几何参数，这对于CT***扫描角度不全的情况非常有针对性。(3)不需要预先知道任何***的几何信息，可以用于自身机械精度很差的几何***，简单方便。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。

在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读介质中，以供指令执行***、装置或设备(如基于计算机的***、包括处理器的***或其他可以从指令执行***、装置或设备取指令并执行指令的***)使用，或结合这些指令执行***、装置或设备而使用。就本说明书而言，"计算机可读介质"可以是任何可以包含、存储、通信、传播或传输程序以供指令执行***、装置或设备或结合这些指令执行***、装置或设备而使用的装置。计算机可读介质的更具体的示例(非穷尽性列表)包括以下：具有一个或多个布线的电连接部(电子装置)，便携式计算机盘盒(磁装置)，随机存取存储器(RAM)，只读存储器(ROM)，可擦除可编辑只读存储器(EPROM或闪速存储器)，光纤装置，以及便携式光盘只读存储器(CDROM)。另外，计算机可读介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为可以例如通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。

应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行***执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。

本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。

此外，在本发明各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。

上述提到的存储介质可以是只读存储器，磁盘或光盘等。尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (11)

1.一种锥束CT成像***的标定方法，其特征在于，包括以下步骤：

建造模体，其中，所述模体为标定参照物，所述模体大体上是塑料圆柱体，所述模体的上下底面各镶嵌了两组钢珠，每组包括N(N不小于5)个钢珠，所述模体底面圆的直径记为2r，所述两组钢珠中心所决定的两个平面的间距记为2l；

以所述钢珠投影在系的所有像素的值为权重对每个钢珠的投影求加权平均坐标值，以获取所述2N个钢珠中心的投影坐标，并分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂；

分别将所述每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序，并得到i系下的原点O和交点D以及计算角度η和A₁、B₁、A₂、B₂四点坐标，其中，所述2N个钢珠中心投影坐标均在所述系下表达，所述A₁、B₁分别为P₁、Q₁在OD上的投影，所述A₂、B₂分别为P₂、Q₂在所述OD上的投影，称矩形P₁P₂Q₁Q₂(长为2l，宽为2r)为探测器绕所述系下的旋转角度θ后，所述模体被i系下的yⁱ-zⁱ平面所截得的截面，所述角度η为所述探测器在绕所述旋转所述角度θ和绕所述系下的旋转角度后，绕所述系下的旋转的角度；

在所述yⁱ-zⁱ平面内，根据所述角度η、所述原点O、所述交点D、Q₁坐标(s₁,t₁)、P₂坐标(s₂,t₂)、A₁坐标A₂坐标B₁坐标B₂坐标构造第一目标函数f(θ*,w_y*,w_z*)，并以此计算所述角度θ、在所述yⁱ-zⁱ平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)，其中，所述原点O_w同时也是所述模体的中心；

根据所述模体两个底面圆上任意一点的投影应分别位于所述轨迹T₁和所述轨迹T₂上的几何关系来构造第二目标函数并计算所述角度

根据所述原点O、所述角度θ、所述角度η、所述角度和在所述yⁱ-zⁱ平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)以及钢珠投影中心坐标间的位置关系构造第三目标函数H(t)并求出角度t；以及

根据所述原点O、所述角度θ、所述角度η、所述角度所述角度t和在所述yⁱ-zⁱ平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)计算在系下的光源P_s坐标和原点O坐标

2.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，

所述系为相对周围环境静止的世界坐标系，所述系以所述模体的中心为原点O_w，所述系以所述模体的轴为在开始投影前，所述系以垂直于所述的射线为所述系下的根据所述所述和右手定则确定；

所述i系为描述所述探测器理想状态的虚拟探测器坐标系，所述i系以所述光源P_s和所述原点O_w的连线与所述探测器平面的交点为所述原点O，所述yⁱ和所述同向，定义P_sO与所述yⁱ确定的平面为所述yⁱ-zⁱ平面，所述zⁱ过所述原点O且垂直于所述yⁱ，xⁱ根据所述yⁱ、所述zⁱ和右手定则确定；以及

所述系为基于所述探测器本身的坐标系，所述系以所述原点O为原点，所述和所述平行于所述探测器的行列线，所述垂直于所述探测器平面，其中，定义所述所述所述分别与所述xⁱ、所述yⁱ、所述zⁱ对应重合时的状态为所述探测器的初始状态。

3.如权利要求2所述的锥束CT成像***的标定方法，其特征在于，从所述探测器的初始状态到所述探测器的任一状态包括以下步骤：

所述探测器从所述探测器的初始状态绕所述旋转，与所述形成所述角度θ，此处，所述角度θ值按照以所述为正方向的右手定则方式定义；

所述探测器从所述探测器的初始状态绕所述旋转，与所述形成所述角度此处，所述角度值按照以所述为正方向的右手定则方式定义；以及

所述探测器从所述探测器的初始状态绕所述旋转，与所述形成所述角度η，此处，所述角度η值按照以所述为正方向的右手定则方式定义。

4.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，根据以下方程分别拟合出每组的N个钢珠中心投影坐标所在椭圆的轨迹T₁和轨迹T₂：

$a_1{(u-u_{0_1})}^2+b_1{(v-v_{0_1})}^2+2c_1(u-u_{0_1})(v-v_{0_1})-1=0$

$a_2{(u-u_{0_2})}^2+b_2{(v-v_{0_2})}^2+2c_2(u-u_{0_2})(v-v_{0_2})-1=0$

其中，a₁、b₁、c₁、和a₂、b₂、c₂、由每组的所述N个钢珠中心投影在所述系的坐标拟合得到。

5.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，所述分别将所述每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序，并得到i系下的原点O和交点D具体为：

分别将所述每组的N个钢珠中心投影坐标在相应的轨迹上定好顺序；

连接E₁F₂和E₂F₁得到所述原点O，其中，所述E₁和所述E₂分别为所述轨迹T₁上纵坐标最大的钢珠中心投影和纵坐标最小的钢珠中心投影，所述F₁和所述F₂分别为所述轨迹T₂上纵坐标最大的钢珠中心投影和纵坐标最小的钢珠中心投影；以及

连接E₁F₁和E₂F₂，所述E₁F₁和所述E₂F₂在空间中延长线的交点作为所述交点D，其中，所述OD与所述的夹角为所述角度η。

6.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，根据以下步骤计算所述角度θ、在所述yⁱ-zⁱ平面内的所述原点O_w坐标(w_y,w_z)和所述光源P_s坐标(s_y,s_z)：

在所述yⁱ-zⁱ平面内，对于待定的原点O_w*坐标(w_y*,w_z*)，根据所述矩形P₁P₂Q₁Q₂的长和宽得到相应的待定矩形顶点P₁*、P₂*、Q₁*、Q₂*，连接所述原点O和所述原点O_w*得到直线OO_w*；

待定角度θ*，过所述交点D作所述yⁱ的平行线DP_s*和所述直线OO_w*相交于点P_s*，连接A₁P₁*、B₁Q₁*、A₂P₂*、B₂Q₂*分别与所述直线OO_w*相交于

构造第一目标函数f(θ*,w_y*,w_z*)用以表示所述P_s*、的离散程度，其中，所述第一目标函数f(θ*,w_y*,w_z*)越大，所述P_s*、的离散程度越大，当所述第一目标函数f(θ*,w_y*,w_z*)＝0时表示所述P_s*、五点重合；

计算 $({\mathrm{\θ}}_0,w_{y_0},w_{z_0})=\arg \min \left\{f(\mathrm{\θ}*,w_y*,w_z*)\right\}$ 的优化结果，以获得所述角度θ和在所述yⁱ-zⁱ平面内的所述原点O_w坐标(w_y,w_z)，即：θ＝θ₀；以及

取对应求出的所述P_s*、五点的平均坐标为所述光源P_s坐标(s_y,s_z)。

7.如权利要求6所述的锥束CT成像***的标定方法，其特征在于，

$f(\mathrm{\θ}*,w_y*,w_z*)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_s*P_{s_i}\left|\right|}^{\mathrm{\λ}},$ 或 $f(\mathrm{\θ}*,w_y*,w_z*)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_s*P_{s_i}\left|\right|}^{\mathrm{\λ}},+\underset{i,j=1,i\≠j}{\overset{4}{\mathrm{\Σ}}}{\left|\right|P_{s_i}{P}_{s_j}\left|\right|}^{\mathrm{\λ}},$ 其中，λ一般为1或2。

8.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，根据以下步骤计算所述角度

待定角度根据已求出的所述角度η、所述角度θ、在所述yⁱ-zⁱ平面内的所述原点O_w坐标(w_y,w_z)、所述光源P_s坐标(s_y,s_z)和待定角度计算所述模体两个底面圆在所述系的轨迹T₁*、T₂*；

在所述轨迹T₁*、T₂*上，分别取点(u₁,v₁)和(u₂,v₂)，构造函数以衡量所述点(u₁,v₁)和(u₂,v₂)分别与所述轨迹T₁和所述轨迹T₂的偏离程度，其中，所述越大，所述点(u₁,v₁)和(u₂,v₂)分别与所述轨迹T₁和轨迹T₂的偏离程度越大；以及

构造第二目标函数以计算所述角度其中，M表示共取了M组在所述轨迹T₁*、T₂*上的所述点(u₁,v₁)和(u₂,v₂)，

9.如权利要求8所述的锥束CT成像***的标定方法，其特征在于，

其中，a₁、b₁、c₁、和a₂、b₂、c₂、由每组的所述N个钢珠中心投影在所述系的坐标拟合得到。

10.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，根据以下公式计算所述角度t：

t＝argmin{H(t*)}

其中，角度参数t*描述所述模体底面圆上的钢珠位于该底面圆上的位置，函数h₁(t*)、h₂(t*)分别为在已知所述角度η、所述角度θ、所述原点O_w坐标(w_y,w_z)，所述光源P_s坐标(s_y,s_z)和所述角度的前提下，在模体两个底面上处于t*所对应位置的钢珠的投影和在所述系下最近的钢珠投影中心的距离。

11.如权利要求1所述的锥束CT成像***的标定方法，其特征在于，