本发明涉及ct技术应用领域，具体说是一种ct系统参数标定及成像算法。
背景技术：
ct是computerizedtomography的缩写，即计算机断层成像。ct技术是指一种在不破坏物体内部结构的前提下，通过从物体周围所获得某种物理量(如波速、x线光强、电子束强等)的投影数据，结合某种数学算法，经过计算机的处理，重建物体特定层面上的二维图像，并根据一系列二维图像构成三维图像的技术。
radon变换是ct技术的主要理论基础。radon变换由数学家radon在1917年提出，radon变换简单来说就是求投影的理论方法。二维情况下，在一个平面内沿与原点距离为d，方向角为θ的直线对原函数f(x,y)做线积分，得到的像函数f(d,θ)就是函数f的radon变换。
随着ct技术的发展与进步，其在物质探测方面所具有的巨大优势日益凸显，使其在医学、工程学、工业、农业、安全检测等领域得到了广泛的应用。但是，ct系统安装时往往存在误差，这影响了成像的质量。因此需要对安装好的ct系统进行参数标定，即借助于已知结构的样品(称为模板)标定ct系统的参数，并据此对未知结构的样品进行成像。
技术实现要素：
本发明的目的是要解决现有技术中存在的不足，提供一种ct系统参数标定及成像算法，用以来对安装好的ct系统进行参数标定。
为达到上述目的，本发明是按照以下技术方案实施的：
一种ct系统参数标定及成像算法，包括解决以下问题：
s1、在正方形托盘上放置两个均匀固体介质组成的椭圆形的标定模板，并在标定模板的两个均匀固体介质上放置两个不同的未知介质的物体，其中标定模板为均匀介质且标定模板的几何信息和吸收率为已知，根据该标定模板及其接收ct系统的数据信息，然后确定ct系统的标定参数：ct系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及x射线的180个方向；
s2、利用二维radon变换和radon反变换对两个未知介质的物体的ct图像重建，确定两个未知介质的物体的几何形状，根据两个未知介质在正方形托盘上的不同位置处的接收ct系统的数据信息来确定两个未知介质的物体的吸收率，进而标定参数及数据确定两个未知介质的物体在正方形托盘中的位置；
具体地，所述s1的具体解决如下：
首先以正方形托盘左下角为原点建立笛卡尔坐标系，对数据进行分析和统计，得出x射线与y轴夹角为0度时的系统状态和x射线与y轴夹角为-90度时的系统状态，通过统计方向为-90度时探测器接收信息的探测单元个数及椭圆介质的几何信息，可求出探测单元之间的距离；
其次，分别统计方向为0度时和方向为-90度时探测器接收信息情况解出旋转中心的位置，分析ct系统的几何结构及该系统绕旋转中心逆时针旋转180次，得出ct系统使用的x射线的180个方向。
具体地，所述s2的具体解决如下：
设f(x)是二维空间的一个函数，radon变换，rf为一个在二维空间中沿着直线的线积分，r为radon变换的算子，公式如下：
rf(l)＝∫lf(x)|dx|
具体的来说，直线l利用长度t来替换，公式：
(x(t),y(t))＝((tsinα+scosα),(-tcosα+ssinα))
其中，s为从原点到直线的l的距离，而α为直线l法向量与x轴之间的角度，因此，
radon变换表达为：
而二维radon变换的几何意义就是函数f(x,y)沿着垂直于直线l的方向的积分，radon反变换则为：
基于傅里叶切片定理的滤波反投影算法：
通过定义物体的二维傅里叶变换来具体阐述傅里叶切片定理的实质内容，如下：
假设一个角度θ的投影为pθ(t)，则它的傅里叶变换为：
当条件θ＝0时，傅里叶切片定理最为简单，首先，当二重积分中v＝0时，考虑是物体在频域上沿着直线的傅里叶变换，则上式的傅里叶变换积分简化为：
根据平行束投影的含义，确定括号里面的项为沿着某个常量x的投影，此时：
由上式可知等式的右边代表投影的一维傅里叶变换，所以，在垂直投影和物体函数的二维变换间满足以下关系：
f(u,0)＝sθ＝0(u)
结合上面的推导，该定理的数学表达式为：
其次，由傅里叶切片定理，推导出基于平行束扫描的重建算法，即滤波反投影算法，傅里叶切片定理表明：一次投影的傅里叶变换是二维傅里叶空间中通过原点的一条直线。
在ct系统中选定固定坐标系，被重建图像为f(x,y)；x射线绕未知介质的物体旋转，因此建立旋转坐标系(s,t)，两坐标系关系为：
t＝xcosθ+ysinθ
s＝-xsinq+ycosq
使用p(t,θ)表示x射线穿过物体衰减的积分，用p(w,θ)表示p(t,θ)的傅里叶变换，根据微积分中的面积元素之间的关系，如下式所示：
其中，j是雅各比行列式，结合傅里叶变换，得到p(w,θ)的表达式
图像二维傅里叶变换表达式为：
令u＝wcosθ，v＝wsinθ，则上面两式相等，并且通过fbp算法整理出f(x,y)表达式，即是平行束扫描重建数学公式，如下式所示：
滤波反投影算法所使用的滤波器为|w|，s-l滤波函数的选取，选用snic函数作为窗函数，得到s-l滤波函数的系统函数：
其相应的冲激响应为：
定义d为投影数据的采样间隔，对应的最高不失真空间频率为b＝1/zd，以t＝1代入上式得到s-l的采样序列的解析表达式为：
用s-l滤波函数重建图像，对该算法及数据作图确定出未知介质的物体的几何形状。
与现有技术相比，本发明具有精度高的特点，得到的结果误差小，接近实际值。
附图说明
图1示出了本发明实施例提供的模板示意图。
图2示出了本发明实施例提供的模板角度分布图。
图3示出了本发明实施例提供的radon变换图。
图4示出了本发明实施例提供的几何意义示意图。
图5示出了本发明实施例提供的radon变换与反变换图。
图6示出了本发明实施例提供的ct系统示意图。
图7示出了本发明实施例提供的fbp算法流程图。
图8示出了本发明实施例提供的未知介质的几何形状示意图。
图9示出了本发明实施例提供的均匀介质衰减示意图。
图10示出了本发明实施例提供的非均匀介质衰减示意图。
图11示出了本发明实施例提供的ct重建图。
具体实施方式
下面结合具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。
为更好地对本发明的算法进行说明，引进一种典型的二维ct系统。在该系统中平行入射的x射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。x射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个x射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后可得到180组接收信息。
第一方面，本发明提供的ct系统参数标定及成像算法，首先解决上述典型的二维ct系统引出的问题一：在正方形托盘上放置两个均匀固体介质组成的标定模板，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。根据这一模板及其接收信息，可确定ct系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该ct系统使用的x射线的180个方向。
第二方面，本发明提供的ct系统参数标定及成像算法，可以解决上述典型的二维ct系统引出的问题二：利用上述典型的二维系统，可以得到某未知介质的接收信息。并可通过已得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。
第三方面，本发明提供的ct系统参数标定及成像算法，可以解决上述典型的二维ct系统引出的问题三：利用上述典型的二维系统，也可以得到另一个未知介质的接收信息。并可通过已得到的标定参数，确定该未知介质的相关信息。
对典型的二维ct系统引出的问题一，进行详细说明：
在正方形托盘左下角建立笛卡尔直角坐标系如图1所示，通过几何图像和接收信息的图像处理求解。
对探测器单元之间的距离的求解：
对从附件二中接收的信息数据进行分析，得到在et列中有156个不为0的数据，即在探测器上有156个探测单元可接收到数据，此为各个方向接收的最少数据，可确定此时为x射线沿着模板椭圆长轴方向时探测器所接收的信息，在bf列有289个不为0的数据，即在探测器上有289个探测单元可接收到数据，此为各个方向接收最多的数据，即可确定此时为bf列为x射线沿着模板椭圆短轴方向时探测器所接收的信息。已知椭圆长轴长为80mm，因此，设探测器单元之间的距离为d，则
故探测器单元之间的距离为0.2778mm。
对ct系统旋转中心的求解：
由附件二接收信息bf列的数据，从探测器92单元到380丹云数据为对应椭圆长轴长接收的信息，所以便有旋转中心y的位置为
由附件二接收信息et列的数据，从探测器169单元到276单元数据为对应椭圆短轴长接收的信息，所以便有旋转中心x的位置为
综上可得，对于此ct系统旋转中心的位置为(40.6944，55.5556)。
对ct系统使用的x射线的180个方向的求解：
通过综合分析附件二et和bf列的数据，et列为x射线沿着椭圆长轴方向时探测器所接收的信息，bf列为x射线沿着模板椭圆短轴方向时探测器所接收的信息，现假设x射线的180个方向用x射线与y轴的夹角表示，整个系统绕旋转中心逆时针旋转，则et列代表的方向为0度，bf列的代表方向为-90度，又在附件二中查得从bf到et列共旋转了93次。从开始到bf列共旋转了57次，故刚开始时的方向为-145.7609度，从et列到结束共旋转了106次，故结束时的方向为29.3478度。综上所述，对此ct系统使用的x射线的180个方向为从-145.7609度到29.3478度。模拟结果如图2所示。
对典型的二维ct系统引出的问题二，进行详细说明：
利用二维radon变换和radon反变换对图像重建，radon反变换运用的是基于傅里叶切片定理的滤波反投影方法，对平行束投影进行重建，确定了未知介质的几何形状，同时得到方针shepp-logan模型，利用改进后算法再实现模型重建。
radon变换：
radon变换后的数据也称为正弦图，而狄拉克函数的radon变换为一个正弦波的分布，因此，物体的radon变换展现为正弦波的组合，如图3所示。
设f(x)是二维空间的一个函数，radon变换，rf为一个在二维空间中沿着直线的线积分，r为radon变换的算子，公式如下：
rf(l)＝∫lf(x)|dx|
具体的来说，直线l可以利用长度t来替换，公式：
(x(t),y(t))＝((tsinα+scosα),(-tcosα+ssinα))
其中，s为从原点到直线的l的距离，而α为直线l法向量与x轴之间的角度。因此，
radon变换可以表达为：
而二维radon变换的几何意义就是函数f(x,y)沿着垂直于直线l的方向的积分，如图4所示。
radon反变换则为：
据此，将radon变换理解为平行束正投影，则其反变换表示平行束ct图像重建过程。做radon变换与反变换结果如图5所示。
基于傅里叶切片定理的滤波反投影算法：
首先，通过定义物体的二维傅里叶变换来具体阐述傅里叶切片定理的实质内容，如下：
假设一个角度θ的投影为pθ(t)，则它的傅里叶变换为：
当条件θ＝0时，傅里叶切片定理最为简单。首先，当二重积分中v＝0时，考虑是物体在频域上沿着直线的傅里叶变换，则上式的傅里叶变换积分可以简化为：
根据平行束投影的含义，可以确定括号里面的项为沿着某个常量x的投影，此时：
由上式可知等式的右边代表投影的一维傅里叶变换，所以，在垂直投影和物体函数的二维变换间满足以下关系：
f(u,0)＝sθ＝0(u)
结合上面的推导，该定理的数学表达式为：
其次，由傅里叶切片定理，我们可以推导出基于平行束扫描的重建算法，即滤波反投影算法。傅里叶切片定理表明：一次投影的傅里叶变换是二维傅里叶空间中通过原点的一条直线。
平行入射的x射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。x射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。ct系统如图6所示。
选定为固定坐标系，因此被重建图像为f(x,y)；x射线绕被测物旋转，因此建立旋转坐标系(s,t)。两坐标系关系为：
t＝xcosθ+ysinθ
s＝-xsinq+ycosq
使用p(t,θ)表示x射线穿过物体衰减的积分，用p(w,θ)表示p(t,θ)的傅里叶变换。根据微积分中的面积元素之间的关系，如下式所示：
其中，j是雅各比行列式。结合傅里叶变换，得到p(w,θ)的表达式
图像二维傅里叶变换表达式为：
令u＝wcosθ，v＝wsinθ，则上面两式相等，并且fbp算法流程图如图7所示。整理出f(x,y)表达式，即是平行束扫描重建数学公式，如下式所示：
滤波反投影算法所使用的滤波器为|w|，s-l滤波函数的选取。为了缓解振荡响应，我们选用snic函数作为窗函数，于是得到s-l滤波函数的系统函数：
其相应的冲激响应为：
定义d为投影数据的采样间隔，对应的最高不失真空间频率为b＝1/zd，以t＝1代入上式得到s-l的采样序列的解析表达式为：
用s-l滤波函数重建图像，因为其振荡相应减小，对含有噪声的投影数据的重建质量较好。对该算法及数据作图确定了未知介质的几何形状如图8所示。
通过对数据进行分析得到几何图形椭圆的长轴长为81.9510mm，短轴长为43.0590mm。该椭圆重心在正方形托盘中的位置是(50，50)且长轴长与y轴成60度夹角。
由于ct重建的数据来源是x射线的衰减，x射线的衰减服从指数规律，均匀介质衰减示意图如图9所示，非均匀介质衰减示意图如图10所示，x射线的入射强度i0、透射强度i和体素的厚度d的关系为：
i＝i0e-μd
因此，物体吸收能量is为：
is＝i0-i＝i0-i0e-μd
然而实际情况是，我们获得的是衰减后的投影数据，因此，投影数据应为：
因为物体为未知介质，吸收率改变，在正方形托盘中位置改变。利用已知标定参数和图像来确定几何形状和在托盘中的位置，进而求解吸收率。未知介质的几何信息结果如图11所示。
当然，为了提高系统参数标定的精度和稳定性，还可以采用归一化平均绝对值距离判断评价标定的参数：
其中，tu，v，γu，v分别表示测试模型和重建后图像中第u行，第v列的吸收能量；
本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。