基于重建图像梯度的L0范数最小化的锥束CT旋转中心标定方法与流程

本发明属于生物医学影像和无损检测技术领域，涉及一种基于重建图像梯度的L0范数最小化的锥束CT(Computed Tomography)旋转中心标定方法。

背景技术：

CT技术成功应用于无损检测领域并取得很大的发展。圆周轨迹的锥束CT系统主要由X射线源，旋转台和探测器三部分组成，锥束CT系统理想的成像关系要求射线源与探测器中心的连线(中心射线)垂直于探测器所在平面，且与旋转台的旋转轴垂直相交，交点就是旋转中心。在对重建图像质量要求高的三维图像重建中，要求锥束CT系统的旋转中心在探测器上的投影位置精确。在实际的锥束CT系统工作中，旋转台和被扫描物体保持不动，射线源和探测器围绕旋转中心一起做圆周运动，由于机械系统的不稳定性，安装完成之后很难保证锥束CT系统实际的旋转中心和理想的旋转中心重合，即理想的旋转中心发生偏移，而在重建算法中默认旋转中心没有发生偏移，与实际情况不符。由于旋转中心的偏移,会导致重建过程中探测器上理想的投影地址和实际的投影地址不符，所以重建图像中会出现几何伪影，使得重建图像的质量不符合精度要求。因此，在图像重建之前需要对锥束CT系统旋转中心的偏移参数进行标定，在重建过程中对旋转中心偏移误差进行补偿，使得重建软件中的旋转中心尽可能接近实际的旋转中心，提高锥束CT系统的成像质量。

为了补偿旋转中心偏移对重建图像质量的影响，Stephen等人于1990年在核科学汇刊上提出了一种利用投影正弦图修正投影中心的方法来计算旋转中心的偏移量，但是这种方法在系统噪声较大时，修正效果并不明显；傅建等人于2003年在兵工学报上提出以焦点的偏差来等效旋转中心的偏移，赋予射线源不同的偏移量，进行图像重建并以信噪比和对比度最高的重建图像对应的偏移量来计算旋转中心真实的偏移量，这种方法选定信噪比和对比度来刻画旋转中心的准确程度，虽然有一定的效果但是标定的精度有限。Tong Liu等人于2006年在光学工程杂志上通过将一根细钢丝固定在旋转台上离旋转中心有一定距离的位置，利用钢丝投影图像的几何关系求出旋转中心的偏移量，该方法将钢丝在探测器上的投影图像理想化为一个点，与实际情况相差较远；李保磊等人于2009年在航空学报上提出一种基于正弦图的工业CT系统旋转中心的标定方法，利用中心投影射线旋转位置不变性，通过计算前后相隔180°的投影数据之差求最小值来确定旋转中心位置；该方法要求射线源剂量非常稳定，否则引入更大误差，使得旋转中心的位置不准确。张蔚等人于2009年在医疗卫生装备杂志上根据投影图关于旋转中心对称原理，先对投影图像进行二值化，再通过求投影图像的质心，最后将多个角度得到的质心坐标求均值得到旋转中心的位置，此方法降低了噪声对求质心位置的影响，但是二值化和求质心产生的误差会严重影响旋转中心位置的确定。

技术实现要素：

有鉴于此，本发明的目的在于提供一种基于重建图像梯度的L0范数最小化的锥束CT旋转中心标定方法，该方法使用的模体制作简单，不需要调整旋转台位置，仅需要一次锥束CT扫描，通过循环微调重建软件中旋转中心的偏移参数，就能够得到较精确的锥束CT系统旋转中心的偏移参数，提高锥束CT系统的成像质量。

为达到上述目的，本发明提供如下技术方案：

一种基于重建图像梯度的L0范数最小化的锥束CT旋转中心标定方法，该方法包括以下步骤：

S1：将均匀的金属圆球放置在旋转台上，射线源和探测器设置在围绕旋转中心的圆形轨道上，启动锥束CT系统，扫描得到金属圆球的投影数据；

S2：根据金属圆球的投影数据，按照锥束CT系统初始的旋转中心位置，通过FDK(Feldkamp，Davis，and Kress)重建算法和重建软件得到金属圆球的重建图像，并计算重建图像梯度的L0范数；

S3：调整重建软件中旋转中心的偏移参数，再次重建得到更新后的重建图像并计算重建图像梯度的L0范数；

S4：循环调整步骤S3中的偏移参数，直到重建图像梯度的L0范数在一定误差范围内达到最小，此时旋转中心的偏移参数就是所要标定的。

进一步，在步骤S1中，锥束CT系统启动时，以旋转中心O为原点，以射线源和旋转中心的连线为Y轴，指向旋转中心的方向为正方向，Z轴垂直于Y轴，并以竖直向上为Z轴正方向，X轴、Y轴与Z轴构成固定的右手笛卡尔坐标系O-XYZ；锥束CT系统启动之后，旋转台在整个扫描过程中保持固定，射线源和探测器一起围绕旋转中心做圆周运动；运动过程中，以旋转中心为原点，以射线源和旋转中心的连线为η轴，指向旋转中心的方向为正方向，Z轴垂直于η轴，ξ轴、η轴与Z轴构成右手笛卡尔坐标系O-ηξZ；X轴与ξ轴的夹角为旋转角θ(0≤θ＜2π)；以探测器中心O_D为原点，U轴和X轴平行且方向一致，V轴和Z轴平行且方向一致，U轴和V轴构成探测器上的右手笛卡尔坐标系O_D-UV；扫描完毕得到金属球的投影数据p_θ(u,v)，p_θ(u,v)表示旋转角度为θ(0≤θ＜2π)时O_D-UV坐标系下坐标为(u,v)的探测器单元上的原始投影数据。

进一步，在步骤S2中，根据金属圆球的投影数据，按照锥束CT系统初始的旋转中心位置，通过FDK重建算法和重建软件得到金属圆球的重建图像，并计算重建图像梯度的L0范数，具体包括以下步骤：

S21：在重建软件中对旋转中心偏移参数进行初始化；

S22：利用FDK重建算法和金属圆球的投影数据重建得到金属圆球的重建图像；

S23：计算步骤S22中重建图像梯度的L0范数。

进一步，在步骤S21中，在重建软件中对旋转中心偏移参数进行初始化，具体包括：锥束CT系统中，旋转中心在Y轴方向上的偏移量只影响重建图像的放大倍数，通常旋转中心距离射线源的距离SO和旋转中心距离探测器的距离SO_D都远大于旋转中心在Y轴方向上的偏移量，因此旋转中心在Y轴方向上的偏移量对放大倍数的影响较小，而旋转中心在X轴和Z轴方向上的偏移量(分别记为横向偏移量Δx和纵向偏移量Δz)对重建图像的影响较大，因此可忽略旋转中心在Y轴方向上的偏移量，重点标定旋转中心的横向和纵向偏移量Δx和Δz；更进一步，旋转中心的横向和纵向偏移量可以等价地转换为探测器中心在U轴和V轴上的偏移量(分别记为横向和纵向偏移量Δu和Δv)，它们之间的关系如下：

因此在利用FDK重建算法时，把初始化旋转中心偏移参数等价为初始化探测器中心的横向和纵向偏移量；定义：

其中，widofDet表示探测器的实际宽度，heiofDet表示探测器的实际高度，ucenter表示横向探测器单元的中心，vcenter表示纵向探测器单元的中心；在重建软件中初始化：

ucenter＝ucenter+Δu₀，vcenter＝vcenter+Δv₀ (3)

其中,Δu₀和Δv₀表示探测器中心的横向和纵向偏移量Δu和Δv的初始值，显然Δu₀和Δv₀可取正值或负值，取正值表示实际的探测器中心沿U轴和V轴正方向偏移，取负值表示实际的探测器中心沿U轴和V轴负方向偏移。

进一步，在步骤S22中，利用FDK重建算法和金属圆球的投影数据重建得到金属圆球的重建图像包括以下几个步骤：

S221：对金属圆球的投影数据进行余弦校正和滤波：

其中，SO表示射线源到旋转中心的距离，u和v分别表示当前探测器单元在坐标系O_D-UV下的横坐标和纵坐标，p_θ(u,v)表示旋转角度为θ时坐标系O_D-UV下坐标为(u,v)的探测器单元上的原始投影数据，表示旋转角度为θ时坐标系O_D-UV下坐标为(u,v)的探测器单元上经过余弦校正和滤波后的投影数据，符号*表示卷积，g(u)表示斜坡滤波器，其在频域的定义为G(w)＝|w|，w表示频域变量；

S222：对滤波后的投影数据进行反投影得到金属圆球的重建图像：

其中，f^FDK(x,y,z)表示重建图像f^FDK在坐标系O-XYZ下坐标为(x,y,z)处的体素值，U(x,y,θ)＝SO-xsin(θ)+ycos(θ)，u(x,y,z,θ)和v(x,y,z,θ)分别表示当前探测器单元在旋转角度为θ时在坐标系O_D-UV下的横坐标和纵坐标，它们可以由体素坐标(x,y,z)和旋转角度θ确定，首先，O-XYZ坐标系下的坐标和O-ηξZ坐标系下的坐标有以下关系：

η＝xcos(θ)+ysin(θ),ξ＝-xsin(θ)+ycos(θ)，z＝z (6)

则有

成立。

进一步，在步骤S23中，计算步骤S22中重建图像梯度的L0范数：在锥束CT图像重建中，重建得到的三维图像可以表示为：

其中，I,J和K表示重建图像在X，Y和Z方向上的维数，i＝0,1,L,I-1，j＝0,1,L,J-1，k＝0,1,L,K-1，表示重建图像在第k层，第j行，第i列处的体素值；重建图像f^FDK在每一点的梯度定义如下：

其中，0＜i＜I-1，0＜j＜J-1,0＜k＜K-1；不失一般性，重建图像外表面体素点的梯度统一设置为0；得到重建图像f^FDK的梯度图像后，统计梯度图像中体素值大于设定阈值T的个数，得到初始状态下重建图像f^FDK梯度的L0范数，记为

进一步，在步骤S3中，调整重建软件中旋转中心的偏移参数，再次重建得到更新后的重建图像并计算重建图像梯度的L0范数，包括以下步骤：

S31：利用(3)式调整重建软件中旋转中心的偏移参数；旋转中心的横向和纵向偏移量可以等价地转换为探测器中心的横向和纵向偏移量；首先，只调整探测器中心的横向偏移量参数Δu，假定其为正值，并赋予Δu一个较大的数值Δu₁；

S32：调用步骤S22和步骤S23得到更新后重建图像梯度的L0范数，记为由于凭借经验选取的Δu₁相对于实际的横向偏移量足够大，因此能够保证不等式成立。

进一步，在步骤S4中，循环调整步骤S3中的偏移参数，直到重建图像梯度的L0范数在一定误差范围内达到最小，此时旋转中心的偏移参数就是所要标定的，具体包括：

S41：在步骤S3的基础上，已知不等式成立，则赋值Δu＝Δu_1/2，其中调用步骤S22和步骤S23得到更新后重建图像梯度的L0范数，记为如果仍有则继续赋值Δu＝Δu_1/4，其中调用步骤S22和步骤S23得到更新后重建图像梯度的L0范数，记为假设时，更新后重建图像梯度的L0范数满足关系其中此时赋值Δu＝Δu_inter，其中，

其中n≤N，N表示对Δu₁逐次减半的最大次数，n表示对Δu₁逐次减半的次数，Δu_inter表示通过线性插值得到的新的Δu值，近似为实际的探测器中心的横向偏移量，记最后的Δu＝Δu_inter，最后利用(1)式将探测器中心的横向偏移量等价地转化为旋转中心的横向偏移量；

S42：在步骤S3的基础上，对Δu₁逐次减半N次数后，其中调用步骤S22和步骤S23得到更新后重建图像梯度的L0范数如果仍有以下不等式成立：

此时，调整重建软件中Δu为负值；

S43：在步骤S42的基础上已知Δu为负值，此时在Δu₁中乘以(-1)，为了保持符号的一致性，仍使用Δu₁，赋值Δu＝Δu₁，并重复步骤S41得到最后的Δu＝Δu_inter，可以近似为实际的探测器中心的横向偏移量，最后利用(1)式将探测器中心的横向偏移量等价地转化为旋转中心的横向偏移量。

本发明的有益效果在于：在本发明中，重建图像梯度的L0范数能够直接反映由旋转中心偏移量导致的几何伪影，并且该方法使用的模体制作简单，锥束CT系统工作中不需要调整旋转台位置，仅需要一次锥束CT扫描，通过循环微调重建软件中旋转中心的偏移参数，就能够得到锥束CT系统较精确的旋转中心偏移量，提高锥束CT系统的成像质量。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明所述方法的流程图；

图2为圆周轨迹的锥束CT系统的示意图；

图3为锥束CT系统工作中，金属圆球上旋转坐标系和固定坐标系之间的关系示意图；

图4为理想的探测器中心和实际的探测器中心之间的关系示意图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

本发明提供的基于重建图像梯度的L0范数最小化的旋转中心偏移参数标定方法，如图1所示，具体包括以下步骤：

步骤1)将均匀的金属圆球5放置在旋转台上，射线源1和探测器3在以旋转中心2为圆心的圆形轨道上，如图2所示，启动锥束CT系统，扫描得到金属圆球的投影数据；锥束CT系统启动时，以旋转中心O为原点，以射线源和旋转中心的连线为Y轴，指向旋转中心的方向为正方向，Z轴垂直于Y轴，并以竖直向上为Z轴正方向，X轴、Y轴与Z轴构成固定的右手笛卡尔坐标系O-XYZ；锥束CT系统启动之后，旋转台在整个扫描过程中保持固定，射线源和探测器一起围绕旋转中心做圆周运动；运动过程中，以旋转中心为原点，以射线源和旋转中心的连线为η轴，指向旋转中心的方向为正方向，Z轴垂直于η轴，ξ轴、η轴与Z轴构成右手笛卡尔坐标系O-ηξZ；X轴与ξ轴的夹角为旋转角θ(0≤θ＜2π)；以探测器中心O_D为原点，U轴和X轴平行且方向一致，V轴和Z轴平行且方向一致，U轴和V轴构成探测器上的右手笛卡尔坐标系O_D-UV；扫描完毕得到金属球的投影数据p_θ(u,v)，p_θ(u,v)表示旋转角度为θ(0≤θ＜2π)时O_D-UV坐标系下坐标为(u,v)的探测器单元上的原始投影数据。

步骤2)根据金属圆球的投影数据，按照锥束CT系统初始的旋转中心位置，通过FDK重建算法和重建软件得到金属圆球的重建图像，并计算重建图像梯度的L0范数；包括以下步骤：

步骤2-1)在重建软件中对旋转中心偏移参数进行初始化；

如图2所示，由于在锥束CT系统中，旋转中心在Y轴方向上的偏移量只影响重建图像的放大倍数，通常旋转中心距离射线源的距离SO和旋转中心距离探测器之间的距离SO_D都远大于旋转中心在Y轴方向上的偏移量，因此旋转中心在Y轴方向上的偏移量对放大倍数的影响较小，而旋转中心在X轴和Z轴方向上的偏移量(分别记为横向偏移量Δx和纵向偏移量Δz)对重建图像的影响较大，因此可忽略旋转中心在Y轴方向上的偏移量，重点标定旋转中心的横向和纵向偏移量Δx和Δz；更进一步，旋转中心的横向和纵向偏移量可以等价地转换为探测器中心在U轴和V轴方向上的偏移量(分别记为横向和纵向偏移量Δu和Δv)，它们之间的关系如下：

因此在利用FDK重建算法时，可以把初始化旋转中心偏移参数等价为初始化探测器中心的横向和纵向偏移量；定义：

其中，widofDet表示探测器的实际宽度，heiofDet表示探测器的实际高度，ucenter表示横向探测器单元的中心，vcenter表示纵向探测器单元的中心；在重建软件中初始化：

ucenter＝ucenter+Δu₀，vcenter＝vcenter+Δv₀ (3)

其中,Δu₀和Δv₀表示探测器中心的横向和纵向偏移量Δu和Δv的初始值，显然Δu₀和Δv₀可取正值或负值，取正值表示实际的探测器中心沿U轴和V轴正方向偏移，取负值表示实际的探测器中心沿U轴和V轴负方向偏移，如图4所示，O'_D-U'V'表示实际的探测器坐标系，O_D-UV表示理想的探测器坐标系，Δu表示实际的探测器中心O'_D相对理想的探测器中心O_D的横向偏移量，Δv表示实际的探测器中心O'_D相对理想的探测器中心O_D的纵向偏移量；

步骤2-2)利用FDK重建算法和金属圆球的投影数据重建得到金属圆球的重建图像；

步骤2-2-1)对金属圆球的投影数据进行余弦校正和滤波：

其中，SO表示射线源到旋转中心的距离，u和v分别表示当前探测器单元在笛卡尔坐标系O_D-UV下的横坐标和纵坐标，p_θ(u,v)表示旋转角度为θ时坐标为(u,v)的探测器单元上的原始投影数据，表示旋转角度为θ时坐标为(u,v)的探测器单元上经过余弦校正和滤波后的投影数据，符号*表示卷积运算，g(u)表示斜坡滤波器，其在频域的定义为G(w)＝w，w表示频域变量；

步骤2-2-2)对滤波后的投影数据进行反投影得到金属圆球的重建图像：

其中，f^FDK(x,y,z)表示重建图像f^FDK在坐标系O-XYZ下坐标为(x,y,z)处的体素值，U(x,y,θ)＝SO-xsin(θ)+ycos(θ)，u(x,y,z,θ)和v(x,y,z,θ)分别表示当前探测器单元在旋转角度为θ时在坐标系O_D-UV下的横坐标和纵坐标，它们可以由物体上体素的坐标(x,y,z)和旋转角度θ确定，首先，O-XYZ坐标系下的坐标和O-ηξZ坐标系下的坐标有以下关系，如图3所示，

η＝xcos(θ)+ysin(θ),ξ＝-xsin(θ)+ycos(θ)，z＝z (6)

则有：

成立。

步骤2-3)计算步骤2-2)中重建图像梯度的L0范数，在锥束CT图像重建中，重建得到的三维图像可以表示为：

其中，I,J和K表示重建图像在X，Y和Z方向上的维数，i＝0,1,L,I-1，j＝0,1,L,J-1，k＝0,1,L,K-1，表示重建图像在第k层，第j行，第i列处的体素值；重建图像f^FDK在每一点的梯度定义如下：

其中，0＜i＜I-1，0＜j＜J-1,0＜k＜K-1；不失一般性，重建图像外表面体素点的梯度统一设置为0；得到重建图像f^FDK的梯度图像后，统计梯度图像中体素值大于设定阈值T的个数，得到初始状态下重建图像f^FDK梯度的L0范数，记为

步骤3)调整重建软件中旋转中心的偏移参数，再次重建得到更新后的重建图像并计算重建图像梯度的L0范数；

步骤3-1)利用(3)式调整重建软件中旋转中心的偏移参数。旋转中心的横向和纵向偏移量可以等价地转换为探测器中心的横向和纵向偏移量。首先，只调整参数Δu，假定其为正值，并赋予Δu一个较大的数值Δu₁；

步骤3-2)调用步骤2-2)和步骤2-3)得到更新后重建图像梯度的L0范数，记为由于凭借经验选取的Δu₁相对于实际的横向偏移量足够大，因此能够保证不等式成立。

步骤4)循环调整旋转中心偏移参数，直到重建图像梯度的L0范数在一定误差范围内达到最小，此时旋转中心的偏移参数就是所要标定的；

步骤4-1)在步骤3)的基础上，已知不等式成立，则赋值Δu＝Δu_1/2，其中调用步骤2-2)和步骤2-3)得到更新后重建图像梯度的L0范数，记为如果仍有则继续赋值Δu＝Δu_1/4，其中调用步骤2-2)和步骤2-3)得到更新后重建图像梯度的L0范数，记为不失一般性，假设时，更新后重建图像梯度的L0范数满足关系其中此时赋值Δu＝Δu_inter，其中，

其中n≤N，N表示对Δu₁逐次减半的最大次数，n表示对Δu₁逐次减半的次数，Δu_inter表示通过线性插值得到的新的Δu值，可以近似为实际的探测器中心的横向偏移量，记最后的Δu＝Δu_inter，最后利用(1)式将探测器中心的横向偏移量等价地转化为旋转中心的横向偏移量；

步骤4-2)在步骤3)的基础上，对Δu₁逐次减半N次数后，其中调用步骤2-2)和步骤2-3)得到更新后重建图像梯度的L0范数如果仍有以下不等式成立：

此时，调整重建软件中Δu为负值；

步骤4-3)在步骤4-2)的基础上已知Δu为负值，此时在Δu₁中乘以(-1)，为了保持符号的一致性，仍使用Δu₁，赋值Δu＝Δu₁，并重复步骤4-1)得到最后的Δu＝Δu_inter，可以近似为实际的探测器中心的横向偏移量，最后利用(1)式将探测器中心的横向偏移量等价地转化为旋转中心的横向偏移量。

需要说明的是，步骤3)和步骤4)首先只考虑了探测器中心的横向偏移量Δu，事实上，计算探测器中心的纵向偏移量Δv和计算探测器中心的横向偏移量Δu有相同的过程，只需要在步骤3)和步骤4)中用符号“v”替换所有的符号“u”，重复步骤3)和步骤4)即可得到探测器中心的纵向偏移量Δv，可以近似为实际的探测器中心的纵向偏移量，最后利用(1)式将探测器中心的纵向偏移量等价地转化为旋转中心的纵向偏移量。

本发明基于重建图像的L0范数最小化去标定锥束CT旋转中心的偏移量，该方法使用的模体制作简单，不需要调整旋转台位置，仅需要一次锥束CT扫描，通过循环微调重建软件中旋转中心的偏移参数，就能够得到较精确的锥束CT系统中旋转中心的偏移量，并重建出清晰的图像。该方法的关键在于使用循环调整旋转中心偏移参数的方法去适应真实的旋转中心位置，而重建图像梯度的L0范数能够直观地反映由旋转中心偏移导致的几何伪影的程度，当重建图像梯度的L0范数在一定误差范围内达到最小时，发明人认为，此时的旋转中心位置和真实的旋转中心位置几乎重合。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其做出各种各样的改变，而不偏离本发明权利要求书所限定的范围。