CN109682843B - 一种对ct系统的参数标定方法 - Google Patents

一种对ct系统的参数标定方法 Download PDF

Info

Publication number
CN109682843B
CN109682843B CN201910116307.0A CN201910116307A CN109682843B CN 109682843 B CN109682843 B CN 109682843B CN 201910116307 A CN201910116307 A CN 201910116307A CN 109682843 B CN109682843 B CN 109682843B
Authority
CN
China
Prior art keywords
data
ray
detector
sample
obtaining
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
Application number
CN201910116307.0A
Other languages
English (en)
Other versions
CN109682843A (zh
Inventor
张雷
向荣艳
张波
张黎
杨国艺
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing Jiaotong University
Original Assignee
Chongqing Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Chongqing Jiaotong University filed Critical Chongqing Jiaotong University
Priority to CN201910116307.0A priority Critical patent/CN109682843B/zh
Publication of CN109682843A publication Critical patent/CN109682843A/zh
Application granted granted Critical
Publication of CN109682843B publication Critical patent/CN109682843B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pulmonology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明公开了一种对CT系统参数进行标定的方法,利用CT系统成像的原理以及射线的衰减规律进行建模,通过已知的模板样品图和180组接收信息数据,对才CT系统的参数进行标定,利用求解的标定参数以某未知介质的接收信息,确定未知介质在托盘中的位置、几何形状,并对标定的结果进行稳定性和精确性分析。本发明建立的模型运算得出相对误差值很小,确保了模型和解决方法的正确性;不仅对CT平行投影的快速重建有着极其便捷的求解思路,而且对各类依据平行光照投影探物有着优良的的适应性,例如精密仪器的光学矫正,各类依靠平行投影的物探,精加工等实物的快速重建都有着重要的意义。也可用于宝石鉴定,辅助传统的影像成像,快速实现鉴定与成像工作。

Description

一种对CT系统的参数标定方法
技术领域
本发明属于CT系统参数标定领域,具体来说,是用已知结构形状的样品(模板)来标定CT系统的各项参数,并对未知样品进行扫描成像的方法。
背景技术
利用物质的吸收特性,在不破坏原物体的情况之下,对生物及工程材料进行断层成像并得到材料的内部信息,这就是CT成像的原理。二维CT成像技术包括平行束入射和锥形束入射。本发明针对的是最简单的平行束入射,其原理如附图1所示,射线平行射向探测器表面,探测器有512个等距的探测器单元,每个单元有一个接收点。射线的发射器和探测器的相对位置固定不变,且从发射器射出X射线与接收器垂直,整个系统绕着固定的旋转中心逆时针旋转180次。对每一个射线的方向,在具有512个等距的单元探测器上测量经位置固定不动的二维待检介质吸收衰减后的射线能量,并经过增益等处理后得到180组接收信息。
在实际情况中,CT系统存在一些误差。安装过程中出现的误差以及CT系统各部件之间的误差,这些误差会影响各项参数和指标的准确性,导致不能与理想系统得出的结果匹配,从而影响成像质量。当误差对系统正常工作影响显著时,则需要对系统进行调整。因此对CT系统进行参数标定尤为重要,用已知结构形状的样品(模板)来标定CT系统的各项参数,并对未知样品进行扫描成像。
发明内容
本发明目的是旨在提供了一种用已知结构形状的样品(模板)来标定CT系统的各项参数,并对未知样品进行扫描成像的方法。
为实现上述技术目的,本发明采用的技术方案如下:
一种对CT系统参数进行标定的方法,包括对探测器单元之间的距离、X射线系统的180个方向、以及旋转中心位置的标定;
其中,对探测器单元之间的距离标定,通过分析模板示意图中的样品形状大小和已知的180组接受信息数据,找出接受到信息数据最多的一组探测器单元,通过测量模板示意图中的信息数据个数结合样品的几何特征以及形状大小来求出探测器单元之间的距离;
对射线系统的180个方向的标定,求出模板示意图中椭圆与圆相切的四条线的切线斜率,计算出切线与坐标的夹角,根据任意两条切线之间的角度差以及其间数据的组数,得出平行射线束的旋转规律,按照旋转规律并结合特殊位置到初始位置之间数据组数,得到初始位置的方向角度,初始位置依次旋转180次得到的位置即为射线的180个方向;
对旋转中心位置的标定,通过射线的衰减规律,可以得出数据信息与切割的椭圆弦长的关系,由于探测器中心射线一定过旋转中心,求出中心射线在不同位置的射线方程,不同位置的方程联立求出的交点就是需要的旋转中心;
对探测器单元之间的距离标定包括如下步骤:
(1)分析模板示意图中的样品形状大小得出椭圆样品的长轴的长度,当X射线在椭圆样品的正右方,探测器在椭圆样品的正左方时,运用MATLAB根据已知的180组接受信息数据所画出的图像,将所得图像与椭圆样品进行比对,数据所成的图像的长度对应着就是椭圆样品的长轴的长度l1
(2)运用MATLAB将方向不同的180组数据中每组最后一次出现接受信息不为0的探测器编号与第一次出现时的探测器编号作差,得到的数值n即方向为该组时穿过椭圆样品时的探测器个数,找出接受到信息数据最多的一组探测器单元,记nmax
(3)将l1和nmax带入公式
d1=l1/nmax (1)
得到探测器单元之间的距离d1
对探测器单元之间的距离标定还包括如下步骤:
(4)当X射线在椭圆样品的正上方,探测器在椭圆样品的正下方时,运用MATLAB根据180组接受信息数据所画出的图像,将所得图像与椭圆样品进行比对,数据所成的图像的长度对应着就是椭圆样品的短轴的长度l2
(5)找出180组数据中接受到信息数据最少的一组探测器单元,记nmin
(6)将l2和nmin带入公式
d2=l2/nmin
得到探测器单元之间的距离d2
(7)将d1和d2进行误差分析,两者相对误差小于1%时,取平均数即得探测器单元之间的距离;
对射线系统的180个方向的标定,包括如下步骤:
(1)通过模板示意图中的样品几何关系,以椭圆中心为原点,建立坐标轴,从而得到模板示意图中椭圆和圆的方程:
Figure GDA0003078155930000031
(x-45)2+y2=42 (3)
(2)假设直线为椭圆和圆的公切线,设直线方程为:
y=kx+m (4)
由于直线与圆相切,可得出方程:
Figure GDA0003078155930000032
(3)根据上述方程可求出四条切线的斜率,之后通过斜率算出切线与坐标轴的角度,在制图软件中画出样品与切线的位置关系,计算出切线与坐标的夹角,根据任意两条切线之间的角度差以及其间数据的组数,得出平行射线束的旋转规律,按照旋转规律并结合特殊位置到初始位置之间数据组数,得到初始位置的方向角度,初始位置依次旋转180次得到的位置即为射线的180个方向;
对旋转中心位置的标定,包括如下步骤:
(1)样品CT数:
Figure GDA0003078155930000033
其中,μ1表示样品的衰减系数,μw表示水的衰减系数;
(2)根据LambertBeer-定律有:
I=I0e-μΔx (7)
其中,I为衰减后接收到的射线强度,I0为入射时X射线的强度,μ为摄像的衰减系数,Δx表示射线切割椭圆的弦长;
由式(7)得到:
Figure GDA0003078155930000041
令L=Δx,t=lnI,t0=lnI0,可以得到一个线性关系式:
Figure GDA0003078155930000042
从长短轴分别垂直X射线所在位置得到的两组数据,将其代入(8)可得出μ.由此可得出增益后数据与射线切割的弦长关系式:
L=μH (10)
其中,H表示增益处理后得到的数据信息;
(3)探测器一共有512束光线,此处选择第256个探测器单元发出的光线作为探测器中心射线,在旋转所得的180个位置,每个位置所在的中心射线都会通过系统的旋转中心,联立不同位置得到的中心射线方程会交于一点,这个点就是旋转中心;
设探测器中心射线的方程为:
y=aix+bi (11)
其中,i=0,1,2,…180,ai为直线的斜率,bi为直线的截距;
将式(11)与上文椭圆公式(2)联立,可得出一个关于x的一元二次方程:
Figure GDA0003078155930000043
将不同位置下的数据带入可求出ai,由韦达定理:
Figure GDA0003078155930000044
在判定公式Δ=b2-4ac≥0的条件下,有弦长公式:
Figure GDA0003078155930000045
(4)根据以上分析,建立几何关系模型如下:
Figure GDA0003078155930000046
代入已知的180组数据并对模型进行求解,可得一系列
Figure GDA0003078155930000051
其中
Figure GDA0003078155930000052
表示直线si的斜率,
Figure GDA0003078155930000053
表示直线si的截距;
选取两条不同位置的中心射线方程,联立求交点即可解出旋转中心
Figure GDA0003078155930000054
作为本发明一种对CT系统参数进行标定的方法的另一种优选,还包括对所标定的参数进行准确性和稳定的分析,对模板示意图中的样品进行几何关系分析,求出特征值的方差,用方差的大小来衡量模型的稳定性。
作为本发明一种对CT系统参数进行标定的方法的另一种优选所述对所标定的参数进行准确性和稳定的分析包括如下步骤:
(1)找出对应椭圆长轴的数据组,将数据代入公式(10)得出一系列的理论弦长值;
(2)利用标定得到的三个标定参数,利用旋转对称的特性,设在不同角度下,通过旋转中心的某条射线所在直线方程为:
y=kix+bi (16)
另设:旋转中心坐标为(x0,y0),则得到:
y0=kix0+b0 (17)
从而解得一系列直线方程,将直线方程与椭圆方程联立,利用前述解得每条直线与椭圆相交所得弦长,再利用其他参数结果,求得弦长的理论值,并对两组结果进行差值方差分析,以此判定模型参数标定的精确性和稳定性。
本发明建立的模型运算得出相对误差值很小,确保了模型和解决方法的正确性,得出的结论也合乎情理;不仅对CT平行投影的快速重建有着极其便捷的求解思路,而且对各类依据平行光照投影探物有着优良的适应性,例如精密仪器的光学矫正,各类依靠平行投影的物探,精加工等实物的快速重建都有着重要的意义。也可用于宝石鉴定,辅助传统的影像成像,快速实现鉴定与成像工作。
附图说明
本发明可以通过附图给出的非限定性实施例进一步说明;
图1为CT系统原理图;
图2为模板示意图,其中单位为mm;
图3为模板示意图中的接收信息;
图4为利用该CT系统,得到的某未知介质的接收信息;
图5为四条切线与模板示意图中椭圆、圆的位置关系图;
图6中心射线相交示意图;
图7为将附图3中的数据反投影图像重建后的结果图;
图8对附图4中数据饭投影图像重建后的图;
图9附图4中未知介质与托盘相对位置图;
图10旋转中心与托盘关系图。
具体实施方式
为了使本领域的技术人员可以更好地理解本发明,下面结合附图和实施例对本发明技术方案进一步说明。
利用图2、3所示,已知的模板样品图和180组接收信息数据(由于数据太多,只列出了1-16组数据中不为0的部分),对才CT系统的参数进行标定,利用求解的标定参数以及图4中的某未知介质的接收信息(数据太多,附图中只列出了1-8组不为0的数据),确定未知介质在托盘中的位置、几何形状和吸收率等信息,并对标定的结果进行稳定性和精确性分析。在本实施例中详细介绍了对探测器单元之间的距离、X射线的方向、旋转中心的坐标的标定的方法模型。
(一)对探测器单元之间的距离标定,通过分析模板示意图中的样品形状大小和已知的180组接受信息数据,找出接受到信息数据最多的一组探测器单元,通过测量模板示意图中的信息数据个数结合样品的几何特征以及形状大小来求出探测器单元之间的距离。
其具体方法如下:
分析附图3中的数据,X射线光源放射出X射线,如果X射线没有穿过椭圆样品,则对应的探测器上接受信息为0,如果吸收信息不为0,则说明X射线在照在探测器的过程中穿过了椭圆样品。
(1)当X射线在椭圆样品的正右方,探测器在椭圆样品的正左方时,运用Matlab根据已知的180组接受信息数据所画出的图像,将所得图像与椭圆样品进行比对,数据所成的图像的长度对应着就是椭圆样品的长轴的长度l1=80mm;
(2)运用Matlab将方向不同的180组数据中每组最后一次出现接受信息不为0的探测器编号与第一次出现时的探测器编号作差,得到的数值n即方向为该组时穿过椭圆样品时的探测器个数,找出接受到信息数据最多的一组探测器单元,记nmax
(3)将l1和nmax带入公式
d1=l1/nmax (1)
得到探测器单元之间的距离d1
通过计算得出第61组数据个数最多,即方向为第61组是穿过椭圆样本的长度最大,所记录的探测器个数nmax=289个,又椭圆中长轴最长。从模板示意图中可知椭圆样品的长轴长度l1=80mm,实际探测器单元之间的距离d1=80/289=0.2768mm。
(4)当X射线在椭圆样品的正上方,探测器在椭圆样品的正下方时,运用MATLAB根据180组接受信息数据所画出的图像,将所得图像与椭圆样品进行比对,数据所成的图像的长度对应着就是椭圆样品的短轴的长度l2=30mm;
(5)找出180组数据中接受到信息数据最少的一组探测器单元,记nmin
(6)将l2和nmin带入公式(1),得到探测器单元之间的距离d2
通过程序计算得出第146组数据个数最少,即方向为第146组是穿过椭圆样本的长度最小,所记录的探测器个数为nmin=197个,又椭圆中短轴轴最短,短轴长度l2=30mm,实际间距d2=30/197=0.2774mm。
(7)将d1和d2进行误差分析,其相对误差系数为:
Figure GDA0003078155930000071
相对误差很小,则所求的两组数值都很接近真实值,于是取平均数:
Figure GDA0003078155930000072
最终得到探测器单元之间的距离为0.2771mm。
(二)对射线系统的180个方向的标定,通过求出模板示意图中椭圆与圆相切的四条线的切线斜率,计算出切线与坐标的夹角,根据任意两条切线之间的角度差以及其间数据的组数,得出平行射线束的旋转规律,按照旋转规律并结合特殊位置到初始位置之间数据组数,得到初始位置的方向角度,初始位置依次旋转180次得到的位置即为射线的180个方向。
具体方法如下:
如前所述,接收数据不为0的数据最少的那组数据,射线基本沿椭圆长轴方向入射,而数据不为0的数量最多的一列则为射线沿短轴方向入射所得数据。
(1)通过模板示意图中的样品几何关系,以椭圆为中心为原点,建立坐标轴,从而得到模板示意图中椭圆和圆的方程:
Figure GDA0003078155930000081
(x-45)2+y2=42 (3)
(2)假设直线为椭圆和圆的公切线,设直线方程为:
y=kx+m (4)
由于直线与圆相切,可得出方程:
Figure GDA0003078155930000082
根据以上分析,我们列出求解切线的几何关系数学模型:
Figure GDA0003078155930000083
(3)对所建立的模型进行求解,得到4条切线的斜率为:
Figure GDA0003078155930000084
之后通过斜率算出切线与坐标轴的角度,在制图软件AutoCAD中画出样品与特殊射线位置,如图5所示。
由图5得知,l1与l2角度相差8°,其间共有8组数据,l2与l3相差89°,其间有89组数据,l3与l1相差86°,其间有86组数据。每转动到一个位置,就会有产生一组数据,由此,基本确定每转一次经过的角度为1°。求出l1与第一列数据之间相差的数据列数为12列,由此得到初始射线入射方位为119°方向。
(三)对旋转中心的坐标的标定。
现有的旋转中心的测量矫正方法可以分为两大类:前处理法与后处理法,它们主要以是否需要进行图像重建分类。前处理法是在重建图像之前对待测样品的原始数据进行分析,判断出旋转中心,此类方法包括几何关系法、对称关系法、中心射线法、针模法、投影边界法等;后处理法是根据重建图像的质量来判断投影中心位置是否正确,如薄壁圆筒法、迭代收敛法等。另外的一种分类方法根据是否需要专用校正模型分为校正模型法与投影原始数据法,如上文所述的所说的针模法与薄壁圆筒法属于前者,几何关系法,投影边界法,中心射线法等属于后者。
本实施例采用前处理法,根据中心射线法与几何关系法相结合,通过椭圆弦长及探测器中心射线,与椭圆方程相结合,求出旋转中心。
具体过程如下:
(1)样品CT数(单位:亨):
Figure GDA0003078155930000091
其中,μ1表示样品的衰减系数,μw表示水的衰减系数(一般取1);
(2)根据LambertBeer-定律有:
I=I0e-μΔx (7)
其中,I为衰减后接收到的射线强度,I0为入射时X射线的强度,μ为摄像的衰减系数,Δx表示射线切割椭圆的弦长;
由式(7)得到:
Figure GDA0003078155930000092
令L=Δx,t=lnI,t0=lnI0,可以得到一个线性关系式:
Figure GDA0003078155930000093
从沿长短轴的特殊位置所得的两组数据,将其代入(8)可得出μ=1.77227,由此可得出增益后数据与射线切割的弦长关系式:
L=μH (10)
其中,H表示增益处理后得到的数据信息;
(3)探测器一共有512束光线,此处选择第256个探测器单元发出的光线作为探测器中心射线,在旋转所得的180个位置,每个位置所在的中心射线都会通过系统的旋转中心,联立不同位置得到的中心射线方程会交于一点,这个点就是旋转中心。
设探测器中心射线的方程为:
y=aix+bi,(i=0,1,2,…180) (11)
其中,ai为直线的斜率,bi为直线的截距;
将式(11)与上文椭圆公式(2)联立,可得出一个关于x的一元二次方程:
Figure GDA0003078155930000101
将不同位置下的数据带入可求出ai,由韦达定理:
Figure GDA0003078155930000102
在判定公式Δ=b2-4ac≥0的条件下,有弦长公式:
Figure GDA0003078155930000103
(4)根据以上分析,建立几何关系模型如下:
Figure GDA0003078155930000104
代入已知的180组数据并对模型进行求解,可得一系列
Figure GDA0003078155930000105
考虑到探测器的入射光线是512条,因此第256条并不是精确的中心射线,这可能会造成一定程度上的误差。基于此,本文选择了10组数据,组成5对中心射线方程组,得到了5个中心坐标位置,我们通过求其平均值,以最后的平均值为真正的旋转中心(如表1所示)。
表1各组中心线的参数表
Figure GDA0003078155930000111
最后,将旋转中心坐标平均,得到最终旋转中心坐标为(-9.6743,5.0399)。
(四)利用上述标定的参数,通过附图4中的接收信息,确定位置介质在托盘中的位置、几何形状。
反投影重建算法的原理是:将样品的切断面上某一点的密度值看作这一平面上所有经过它的所有射线投影之和。反投影重建算法是通过寻找经过某个像素点的所有射线投影,将每个射线投影沿射线方向均匀地分配投影值,根据反投影结果映射生成图像。但对于CT成像系统,一个理想点状的物体都会扩散成为一个分布,形成星状伪影,造成重建图像的严重失真。为去除形状伪影,就需要在反投影重建之前对数据进行滤波处理。
平行束滤波反投影重建算法则是根据反投影算法演变而来,它是通过傅里叶切片定理推导得到,待重建的图像为f(xy),则其滤波反投影重建公式为:
Figure GDA0003078155930000112
其中,h(xr)与p(xr,α)都由傅里叶切片定理得到。
傅里叶切片定理:
Figure GDA0003078155930000113
其中,G1[·]表示一维傅里叶变换,G(ρ,α)是二维傅里叶变换的极坐标表示方法。而pα(xr)为f(x,y)在α=α0的平行束投影。
平行线滤波反投影重建图像的具体过程是,先将探测器上得到的数据与滤波器函数进行卷积运算,得到射线不同方向上的投影数据;而后将数据沿每个方向进行反投影(按原来的方向平均分配到每一个矩阵单元,叠加后即可得到每组矩阵单元的CT值),经过适当处理,原物像被还原。
算法步骤如下:
step1:设计合适的滤波器h;
step2:把在角度αi时测到的投影p(xr,αi)与滤波器的h进行卷积滤波,就求出了经过滤波后的投影g(xr,αi);
step3:对每个αi,使其滤波后得到的投影p(xr,αi)的反投影满足xr=rcos(β-α)的射线上的所有点r,β;
step4:将step3中算出的反投影值对所有的0≤αi≤2π进行累加,就能够得到重建后的图像数据。
将附图3中的数据导入Matlab,根据反投影算法与经滤波后的反投影算法程序运行,我们得到结果如图7所示。
通过重建图像,将图7与图2对比,可以知道滤波后的成图效果明显比未经滤波处理的效果好。由此,根据滤波反投影,我们在Matlab中运行出附图4的图形如图8所示。
由图2以及已知的数据信息分析得知,托盘的正方向与入射光线之间存在一个夹角为119°,自左上向右下入射,根据假设,中心光线始终穿过旋转中心可知在任意旋转角度下始终有512条光束的中心光束通过旋转中心,并且旋转中心与托盘间始终存在一定的几何关系,如图9(外框为托盘)。其中,点(x0,y0)由模型一的建立求解可得为(40.33,44.99)。
基于以上分析可依据入射X射线与托盘间的夹角以及旋转中心与托盘间的位置关系对Matlab的求解得到的重建的图形进行位置的确定,具体如下:
(1)将图形导入AutoCAD中调整重建图形的大小,将其边长调整为为141.42mm的正方形水平放置,并连接其两条对角线,焦点即为旋转中心。
(2)将重建图形绕着旋转中心逆时针旋转29度得到托盘的位置范围区间。
(3)根据图9的旋转中心与托盘的位置关系得出托盘的具体位置如图10中,最内侧方框所示位置。
(五)对前述标定的参数的稳定性和精确性的分析。
(1)找出对应椭圆长轴的数据组,将数据代入公式(10):L=1.77227H,得出一系列的理论弦长值;
(2)利用标定得到的三个标定参数,利用旋转对称的特性,设在不同角度下,通过旋转中心的某条射线所在直线方程为:
y=kix+bi (16)
另设:旋转中心坐标为(x0,y0),则得到:
y0=kix0+b0 (17)
从而解得一系列直线方程,将直线方程与椭圆方程联立,利用前述解得每条直线与椭圆相交所得弦长,再利用其他参数结果,求得弦长的理论值,并对两组结果进行差值方差分析,以此判定模型参数标定的精确性和稳定性。
方差公式(这里任意取10组数据):
Figure GDA0003078155930000131
得出如下数据表。
表2弦长差值分析表
Figure GDA0003078155930000132
从表格中数据求得差值的均值为:
Figure GDA0003078155930000133
解出s2=0.8060,比较小,在可以接受的范围内。因此,本文对参数标定所建立的一系列几何关系模型是比较准确的。
需要说明的是,本实施例提供的方法模型是在X射线的能量只有在穿过实物才会造成能量衰减;旋转中心穿过的是第256个探测器单元发出的x射线的假设下成立的。
以上对本发明提供的一种CT系统的参数标定及利用该参数成像的方法进行了详细介绍。具体实施例的说明只是用于帮助理解本发明的方法及其核心思想。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以对本发明进行若干改进和修饰,这些改进和修饰也落入本发明权利要求的保护范围内。

Claims (3)

1.一种对CT系统参数进行标定的方法,其特征在于,包括对探测器单元之间的距离、X射线系统的180个方向、以及旋转中心位置的标定;
其中,对探测器单元之间的距离标定,通过分析模板示意图中的样品形状大小和已知的180组接受信息数据,找出接受到信息数据最多的一组探测器单元,通过测量模板示意图中的信息数据个数结合样品的几何特征以及形状大小来求出探测器单元之间的距离;
对射线系统的180个方向的标定,求出模板示意图中椭圆与圆相切的四条线的切线斜率,计算出切线与坐标的夹角,根据任意两条切线之间的角度差以及其间数据的组数,得出平行射线束的旋转规律,按照旋转规律并结合特殊位置到初始位置之间数据组数,得到初始位置的方向角度,初始位置依次旋转180次得到的位置即为射线的180个方向;
对旋转中心位置的标定,通过射线的衰减规律,可以得出数据信息与切割的椭圆弦长的关系,由于探测器中心射线一定过旋转中心,求出中心射线在不同位置的射线方程,不同位置的方程联立求出的交点就是需要的旋转中心;
对探测器单元之间的距离标定包括如下步骤:
(1)分析模板示意图中的样品形状大小得出椭圆样品的长轴的长度,当X射线在椭圆样品的正右方,探测器在椭圆样品的正左方时,运用MATLAB根据已知的180组接受信息数据所画出的图像,将所得图像与椭圆样品进行比对,数据所成的图像的长度对应着就是椭圆样品的长轴的长度l1
(2)运用MATLAB将方向不同的180组数据中每组最后一次出现接受信息不为0的探测器编号与第一次出现时的探测器编号作差,得到的数值n即方向为该组时穿过椭圆样品时的探测器个数,找出接受到信息数据最多的一组探测器单元,记nmax
(3)将l1和nmax带入公式
d1=l1/nmax (1)
得到探测器单元之间的距离d1
对探测器单元之间的距离标定还包括如下步骤:
(4)当X射线在椭圆样品的正上方,探测器在椭圆样品的正下方时,运用MATLAB根据180组接受信息数据所画出的图像,将所得图像与椭圆样品进行比对,数据所成的图像的长度对应着就是椭圆样品的短轴的长度l2
(5)找出180组数据中接受到信息数据最少的一组探测器单元,记nmin
(6)将l2和nmin带入公式
d2=l2/nmin
得到探测器单元之间的距离d2
(7)将d1和d2进行误差分析,两者相对误差小于1%时,取平均数即得探测器单元之间的距离;
对射线系统的180个方向的标定,包括如下步骤:
(1)通过模板示意图中的样品几何关系,以椭圆中心为原点,建立坐标轴,从而得到模板示意图中椭圆和圆的方程:
Figure FDA0003078155920000021
(x-45)2+y2=42 (3)
(2)假设直线为椭圆和圆的公切线,设直线方程为:
y=kx+m (4)
由于直线与圆相切,可得出方程:
Figure FDA0003078155920000022
(3)根据上述方程可求出四条切线的斜率,之后通过斜率算出切线与坐标轴的角度,在制图软件中画出样品与切线的位置关系,计算出切线与坐标的夹角,根据任意两条切线之间的角度差以及其间数据的组数,得出平行射线束的旋转规律,按照旋转规律并结合特殊位置到初始位置之间数据组数,得到初始位置的方向角度,初始位置依次旋转180次得到的位置即为射线的180个方向;
对旋转中心位置的标定,包括如下步骤:
(1)样品CT数:
Figure FDA0003078155920000023
其中,μ1表示样品的衰减系数,μw表示水的衰减系数;
(2)根据LambertBeer-定律有:
I=I0e-μΔx (7)
其中,I为衰减后接收到的射线强度,I0为入射时X射线的强度,μ为摄像的衰减系数,Δx表示射线切割椭圆的弦长;
由式(7)得到:
Figure FDA0003078155920000024
令L=Δx,t=lnI,t0=lnI0,可以得到一个线性关系式:
Figure FDA0003078155920000031
从长短轴分别垂直X射线所在位置得到的两组数据,将其代入(8)可得出μ,由此可得出增益后数据与射线切割的弦长关系式:
L=μH (10)
其中,H表示增益处理后得到的数据信息;
(3)探测器一共有512束光线,此处选择第256个探测器单元发出的光线作为探测器中心射线,在旋转所得的180个位置,每个位置所在的中心射线都会通过系统的旋转中心,联立不同位置得到的中心射线方程会交于一点,这个点就是旋转中心;
设探测器中心射线的方程为:
y=aix+bi (11)
其中,i=0,1,2,…180,ai为直线的斜率,bi为直线的截距;
将式(11)与上文椭圆公式(2)联立,可得出一个关于x的一元二次方程:
Figure FDA0003078155920000032
将不同位置下的数据带入可求出ai,由韦达定理:
Figure FDA0003078155920000033
在判定公式Δ=b2-4ac≥0的条件下,有弦长公式:
Figure FDA0003078155920000034
(4)根据以上分析,建立几何关系模型如下:
Figure FDA0003078155920000035
代入已知的180组数据并对模型进行求解,可得一系列
Figure FDA0003078155920000036
其中
Figure FDA0003078155920000037
表示直线si的斜率,
Figure FDA0003078155920000038
表示直线si的截距;
选取两条不同位置的中心射线方程,联立求交点即可解出旋转中心
Figure FDA0003078155920000041
2.根据权利要求1所述的一种对CT系统参数进行标定的方法,其特征在于,还包括对所标定的参数进行准确性和稳定的分析,对模板示意图中的样品进行几何关系分析,求出特征值的方差,用方差的大小来衡量模型的稳定性。
3.根据权利要求2所述的一种对CT系统参数进行标定的方法,其特征在于,所述对所标定的参数进行准确性和稳定的分析包括如下步骤:
(1)找出对应椭圆长轴的数据组,将数据代入公式(10)得出一系列的理论弦长值;
(2)利用标定得到的三个标定参数,利用旋转对称的特性,设在不同角度下,通过旋转中心的某条射线所在直线方程为:
y=kix+bi (16)
另设:旋转中心坐标为(x0,y0),则得到:
y0=k0x0 +b0 (17)
从而解得一系列直线方程,将直线方程与椭圆方程联立,利用前述解得每条直线与椭圆相交所得弦长,再利用其他参数结果,求得弦长的理论值,并对两组结果进行差值方差分析,以此判定模型参数标定的精确性和稳定性。
CN201910116307.0A 2019-02-13 2019-02-13 一种对ct系统的参数标定方法 Active CN109682843B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910116307.0A CN109682843B (zh) 2019-02-13 2019-02-13 一种对ct系统的参数标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910116307.0A CN109682843B (zh) 2019-02-13 2019-02-13 一种对ct系统的参数标定方法

Publications (2)

Publication Number Publication Date
CN109682843A CN109682843A (zh) 2019-04-26
CN109682843B true CN109682843B (zh) 2021-07-06

Family

ID=66195872

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910116307.0A Active CN109682843B (zh) 2019-02-13 2019-02-13 一种对ct系统的参数标定方法

Country Status (1)

Country Link
CN (1) CN109682843B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110243847A (zh) * 2019-07-04 2019-09-17 湖南理工学院 一种ct系统参数标定及成像方法

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4464775A (en) * 1980-10-08 1984-08-07 Tokyo Shibaura Denki Kabushiki Kaisha Method and apparatus for collecting X-ray absorption data in X-ray computed tomography apparatus
CN1861001A (zh) * 2006-03-31 2006-11-15 北京航空航天大学 一种适用于三维ct扫描系统投影坐标原点的标定方法
CN101061957A (zh) * 2006-04-28 2007-10-31 西门子公司 在x射线ct中进行散射校正的方法和应用该方法的ct
DE102010042388A1 (de) * 2010-10-13 2012-04-19 Siemens Aktiengesellschaft Verfahren zur Kalibrierung eines CT-Systems mit zumindest einer Fokus-Detektor-Kombination mit einem quantenzählenden Detektor
WO2014154627A1 (en) * 2013-03-27 2014-10-02 Nikon Metrology Nv Registration object, correction method and apparatus for computed radiographic tomography
CN105092612A (zh) * 2014-05-06 2015-11-25 阿斯特菲公司 用于货物和运输容器的计算机断层扫描系统
CN105849537A (zh) * 2013-11-28 2016-08-10 尼康计量公众有限公司 校准设备和计算机断层扫描方法
CN108182665A (zh) * 2017-12-26 2018-06-19 南京邮电大学 一种基于滤波反投影‐迭代算法的ct系统图像重建方法
CN108333197A (zh) * 2018-02-01 2018-07-27 北京航空航天大学 偏置扫描模式下工业ct系统转台旋转中心标定方法
CN109157238A (zh) * 2018-10-19 2019-01-08 中国科学院长春光学精密机械与物理研究所 一种提高x射线ct成像精度的ct数据校正方法
CN109765249A (zh) * 2018-12-26 2019-05-17 陕西师范大学 一种基于random变换的CT系统参数标定与成像方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011143239A (ja) * 2009-12-14 2011-07-28 Toshiba Corp X線ct装置及びその制御方法
US8929507B2 (en) * 2011-10-19 2015-01-06 Kabushiki Kaisha Toshiba Method and system for substantially reducing ring artifact based upon ring statistics

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4464775A (en) * 1980-10-08 1984-08-07 Tokyo Shibaura Denki Kabushiki Kaisha Method and apparatus for collecting X-ray absorption data in X-ray computed tomography apparatus
CN1861001A (zh) * 2006-03-31 2006-11-15 北京航空航天大学 一种适用于三维ct扫描系统投影坐标原点的标定方法
CN101061957A (zh) * 2006-04-28 2007-10-31 西门子公司 在x射线ct中进行散射校正的方法和应用该方法的ct
DE102010042388A1 (de) * 2010-10-13 2012-04-19 Siemens Aktiengesellschaft Verfahren zur Kalibrierung eines CT-Systems mit zumindest einer Fokus-Detektor-Kombination mit einem quantenzählenden Detektor
WO2014154627A1 (en) * 2013-03-27 2014-10-02 Nikon Metrology Nv Registration object, correction method and apparatus for computed radiographic tomography
CN105849537A (zh) * 2013-11-28 2016-08-10 尼康计量公众有限公司 校准设备和计算机断层扫描方法
CN105092612A (zh) * 2014-05-06 2015-11-25 阿斯特菲公司 用于货物和运输容器的计算机断层扫描系统
CN108182665A (zh) * 2017-12-26 2018-06-19 南京邮电大学 一种基于滤波反投影‐迭代算法的ct系统图像重建方法
CN108333197A (zh) * 2018-02-01 2018-07-27 北京航空航天大学 偏置扫描模式下工业ct系统转台旋转中心标定方法
CN109157238A (zh) * 2018-10-19 2019-01-08 中国科学院长春光学精密机械与物理研究所 一种提高x射线ct成像精度的ct数据校正方法
CN109765249A (zh) * 2018-12-26 2019-05-17 陕西师范大学 一种基于random变换的CT系统参数标定与成像方法

Also Published As

Publication number Publication date
CN109682843A (zh) 2019-04-26

Similar Documents

Publication Publication Date Title
Lawrence et al. Transform-based backprojection for volume reconstruction of large format electron microscope tilt series
CN102711613B (zh) 计算断层摄影成像方法及系统
Meng et al. Online geometric calibration of cone-beam computed tomography for arbitrary imaging objects
CN111612768B (zh) 一种采用结构光空间定位和二维工业ct检测叶片方法
EP3814759B1 (en) Item inspection by radiation imaging using an iterative projection-matching approach
CN109949411B (zh) 一种基于三维加权滤波反投影和统计迭代的图像重建方法
CN102652674A (zh) 一种消除ct图像中的几何伪影的方法和系统
US20080212734A1 (en) Correction of Non-Linearities in an Imaging System by Means of a Priori Knowledge in Radiography
CN111553849A (zh) 基于局部特征匹配的锥束ct几何伪影去除方法及装置
CN105510362A (zh) 基于微型ct的水稻分蘖性状无损测量装置及其测量方法
CN103479379A (zh) 一种倾斜螺旋扫描的图像重建方法及装置
CN109682843B (zh) 一种对ct系统的参数标定方法
Yang et al. Geometry calibration method for a cone‐beam CT system
CN107016655A (zh) 锥束cl几何全参数迭代校正方法
CN117874900A (zh) 一种基于bim技术的房屋建筑工程监理方法
Duan et al. Knowledge-based self-calibration method of calibration phantom by and for accurate robot-based CT imaging systems
JP2023005078A (ja) 角度誤差推定装置、方法およびプログラム
CN109884090B (zh) 一种改进圆盘卡法的ct空间分辨率测量方法
EP2453798B1 (en) System and method for image reconstruction by using multi-sheet surface rebinning
CN109765249A (zh) 一种基于random变换的CT系统参数标定与成像方法
Nguyen et al. BeadNet: a network for automated spherical marker detection in radiographs for geometry calibration
US7702180B2 (en) Imaging method and device for the computer-assisted evaluation of computer-tomographic measurements by means of direct iterative reconstruction
CN106910163A (zh) 原始ct投影数据的数据修复装置及方法、ct成像系统
CN110974280A (zh) 一种适用于锥束ct的几何参数精确自动化校正方法
CN118334162B (zh) 激光增材制造板状构件x射线层析成像几何参数标定方法

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
GR01 Patent grant
GR01 Patent grant