一种基于random变换的CT系统参数标定与成像方法与流程

本发明属于ct技术应用领域，具体涉及一种基于random变换的ct系统参数标定与成像方法。

背景技术：

ct技术利用x射线穿透物体的衰减信息重建出物体断层图像。现考虑一个二维ct系统：在探测器平面垂直方向，x射线平行入射，该探测器的512个探测器单元等距排列。固定二维待检测介质的位置，让整个发射-接收系统绕某固定的旋转中心逆时针旋转180次，换言之，即x射线方向改变180次。二维待检测介质吸收衰减后的射线能量可测得，经过增益等处理后得到180组接收信息。为了提高成像质量，往往借助于模板来标定ct系统的参数，并可据此对未知结构的样品进行成像。

技术实现要素：

为了解决现有技术中存在的成像质量欠佳的问题，本发明提供了一种基于random变换的ct系统参数标定与成像方法。本发明要解决的技术问题通过以下技术方案实现：

一种基于random变换的ct系统参数标定与成像方法，包括：

问题1、在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息和吸收率已知，根据该标定模板及其接收信息，求解ct系统旋转中心在正方形托盘中的位置、探测器单元的间距以及x射线的180个方向；

问题2、利用问题1中ct系统得到的未知介质1的接收信息，未知介质1相应的数据文件已知，利用问题1中的标定参数，确定该未知介质1在正方形托盘中的位置、几何形状及吸收率；

问题3、利用问题1中ct系统得到的另一个未知介质2的接收信息，未知介质2相应的数据文件已知，利用问题1中的标定参数，给出该未知介质2在正方形托盘中的位置、几何形状及吸收率；

问题4、设计新模板建立对应的标定模型，以改进ct系统的标定精度和稳定性。

进一步地，上述ct系统的参数标定的具体方法如下：

当被测物体是均匀介质时，穿过物体的射线强度按照指数规律衰减，遵循lambert-beer定律，用公式表示出来就是w＝w0e^-αl；

而在射线只通过空气的路径上，x射线的强度几乎没有改变，因此，空气的值几乎为零；

在x-y平面内，若衰减系数的分布函数为p(x,y)，则在l所在直线的方向上的射线强度变化为

取负对数后，记为f，是一个可测物理量，它的物理意义是物质在该能量射线下的衰减系数沿直线l方向上的线积分，

∫lp(x,y)dl＝αln(w/w0)

由于本题中介质均匀，因此分布函数p(x,y)恒为定值。

进一步地，上述未知介质1的位置、几何形状和吸收率的确定方法如下：

ct成像相当于对一切探测出的投影值rg(t,θ)，反求衰减系数的分布函数g(x,y)，这个反演过程可通过拉东变换，

拉东变换的参数形式定义式为

其中g(x,y)函数沿平面上任何一条直线xcosθ+ysinθ＝t可积，上式称为函数g(x,y)的拉东变换，记b(t,θ)＝rg(t,θ)；

先求出ct系统的增益系数，然后对数据进行去增益化

其中α＝ln[(l0-l)l]，计算出α＝1.75；

由于ct系统的旋转中心不在正方形托盘的中心，因此需要对数据做平移变化才能恢复数据的原始信息；

设旋转中心为(a,b)，做如下平移变化，

rρ(t,θ)＝rf(t-acosθ-bsinθ,θ)

由中心切片定理知，图像重构就是求投影g(t,θ)的逆变换(反演变换)，对当θ给定时，图像投影rg(t,θ)的一维傅里叶变换等于对图像g(x,y)的二维傅里叶变换；

首先对图像投影做一维傅里叶变换，

再做二维傅里叶变换，

此时通过观测值b(t,θ)反推出的f(x,y)为经过校正后的信息。

进一步地，上述未知介质2的位置、几何形状和吸收率的确定方法与未知介质1的相关信息的理论求解过程一致，相互区别的只是相应的经过增益后的数据不同。

进一步地，上述ct系统的误差分析与稳定性分析的具体方法如下：在ct系统的旋转次数是有限的情况下，图像重建是不适定的，利用被测物体的先验信息对问题进行约束，可求出原问题的近似解，正则化算法就是这样一种方法；

正则化方法的求解过程为：在对线性代数方程组ag＝b进行求解时，可将此问题转化为计算||ag-b||2的极小值，但对于不适定问题的求解，可添加一个惩罚项，此时原问题变为：

其中，ψ(g)取为正则函数，为附加在正则函数上的惩罚项因子，对极小值进行求解也就是对公式最优解进行求解；

根据几何关系，可以计算出探测点元之间的距离，

当斜率为定值时，x1/x2也为一定值，由此可以计算出r的系数为(x1+1/x2)，因为r和x1/x2代入算出距离相对较小，由此可以判断小圆半径对稳定性的影响较小。

与现有技术相比，本发明的有益效果：

1.本发明由衰减信息的增减变化与射线穿过物体的投影的对应关系，利用判别式法、几何法、中垂线法建立关于斜率与截距的ct系统旋转角度模型，从优化模型的角度出发提高参数的精度，对数据进行校正处理，利用matlab自带的wdencmp函数进行小波去噪处理，得到ct图像，成像质量高；

2.本发明建立的模型利用几何方法求解出了ct系统的各个参数，直观易懂；同时，在计算ct系统的旋转中心时，随机选取的数据计算出的参数已经达到一定精度，数据量不大，计算相对简单，是一个很好的求近似解的办法。

附图说明

图1是本发明ct系统示意图。

图2是本发明模板示意图。

图3是本发明10个位置示意图。

图4是本发明ct成像的投影示意图。

图5是本发明投影信息不重合示意图。

图6是本发明投影信息重合示意图。

图7是本发明射线为0时的示意图。

图8是本发明投影信息重合示意图。

图9是本发明投影信息例图。

图10是本发明180个相邻角度的差值变化图。

图11是本发明中垂线法示意图。

图12是本发明未知介质1的ct图像。

图13是本发明未知介质2的ct图像。

图14是本发明精度分析例图。

图15是本发明第四象限x射线示意图。

具体实施方式

下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。

为了解决现有技术中存在的成像质量欠佳的问题，本实施例提供了一种基于random变换的ct系统参数标定与成像方法，参照图1-图15，该基于random变换的ct系统参数标定与成像方法，包括：

问题1、在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息和吸收率已知，模板的几何信息如图2，根据该标定模板及其接收信息，求解ct系统旋转中心在正方形托盘中的位置、探测器单元的间距以及x射线的180个方向；

问题2、利用问题1中ct系统得到的未知介质1的接收信息，利用问题1中的标定参数，确定该未知介质1在正方形托盘中的位置、几何形状及吸收率，并给出图3所给的10个位置处的吸收率；

问题3、利用问题1中ct系统得到的另一个未知介质2的接收信息，未知介质2相应的数据文件已知，利用问题1中的标定参数，给出该未知介质2在正方形托盘中的位置、几何形状及吸收率，并给出图3所给的10个位置处的吸收率；

问题4、设计新模板建立对应的标定模型，以改进原ct系统的标定精度和稳定性。

该二维ct系统的x射线只穿过了两种介质，一个是待测介质，另一个是空气。一方面，由于空气引起的衰减很少，射线强度衰减程度小，因此吸收率为0。另一方面，由于该二维待测介质是均匀的，所以介质上各点的吸收率为1。这与模板每一点的吸收率的数值是符合的。射线强度的衰减情况通常与通过物质的厚度、密度、成分有着密切联系，但本题中ct系统与待测介质均是二维的，且待测介质密度均匀，因此只需考虑射线在介质中经过的路径造成的衰减。

针对问题一，首先，考虑以椭圆的中心为原点，以短轴方向为x轴，长轴方向为y轴建立平面直角坐标系，通过射线中的特殊直线，如椭圆切线，联立直线方程与椭圆方程，再寻找图中可能存在的某些数量关系，求解出直线斜率，即可求出180个旋转角度。其次，在求解旋转角度的过程中，已经求出了在射线方向确定时，椭圆的切线方程。当一族平行直线穿过一关于原点对称的椭圆时，经过原点的直线被椭圆所截的线段最长。模板的接收信息是180个方向下，探测器的512个单元所接收到的信息，大致观察这些信息的增减趋势可以了解到，每一列中的最大值可能对应于经过椭圆中心，即原点的射线，而当数据被一段零值阻隔后又重新有非零数值就是代表在这个方向上有射线经过椭圆和圆之间的空气部分。在给定射线方向时，平行直线与椭圆有且仅有两条切线，这两条切线之间的距离也是可计算的，再根据模板的接收信息可以找出两条切线中间占了多少探测器单元，由于探测器单元等距排列，因此切线间的距离与对应探测器单元总数的比值即为所求的探测器单元间的间距。最后，不论ct系统旋转多少角度，射线方向如何，可以肯定的是，旋转中心一定在探测器所在线段的中垂线上，两条不平行直线确定一个交点，因此可以通过在模板的接收信息中任取2个方向的数据，用中垂线法求出旋转中心g的坐标。

针对问题二，在ct系统参数已标定的情况下，如何确定ct系统得到的未知介质1的接收信息对应的未知介质在正方形托盘中的位置、几何形状和吸收率等信息，换言之，就是已知观测值时，在上述ct系统中如何确定模板信息的问题。由于ct系统初始安装时总会存在误差，且实际测量出的接收信息是经过ct系统的增益处理的，因此需要对观测值进行校正，利用校正后的数据通过拉东变换的重建图像，将得到所求的模板信息与接收信息。

针对问题三，与问题二所需求解的内容一样，不同的是待测介质不同，接收到的数据也不同。重复问题二的过程应可求解出该未知介质的相关信息与图3所给的10个位置处的吸收率。

针对问题四，若要提高ct系统成像的精度，大致有三个方面，增大射线的密集程度，模型是否最优，校正ct系统的参数。由于该二维ct系统的探测单元数为定值512，且各探测单元等距排列，因此射线的密集程度是确定的，无法改变的，考虑从对模型进行优化的角度出发，求出模型的最优近似解以减少误差。正则化方法的本质是加入一些先验信息进行约束，得到原问题的近似解。

模型假设：假设x射线是单能的；由空气引起的射线强度的衰减很小；不考虑周围环境的辐射干扰；不考虑待测物质自身产生的辐射；待测物质分别为标准椭圆和标准圆。

一、ct系统的参数标定的具体方法如下：

当被测物体是均匀介质时，穿过物体的射线强度按照指数规律衰减，遵循lambert-beer定律，用公式表示出来就是w＝w0e^-αl；

而在射线只通过空气的路径上，x射线的强度几乎没有改变，因此，空气的值几乎为零；

在x-y平面内，若衰减系数的分布函数为p(x,y)，则在l所在直线的方向上的射线强度变化为

取负对数后，记为f，是一个可测物理量，它的物理意义是物质在该能量射线下的衰减系数沿直线l方向上的线积分，

∫lp(x,y)dl＝αln(w/w0)

由于本题中介质均匀，因此分布函数p(x,y)恒为定值。观察已知数据只有0,1两种数值，因此该定值是1。

1)ct系统的旋转角度模型一

取定x射线的方向如图1所示，按图示方向考虑接收器的接受信息的变化，随着l的增加，衰减信息由零增加到峰值，再减少，若存在射线穿过椭圆与圆的空隙，即射线仅穿过空气，则衰减信息在之前的基础上减少到零，并可能在一段区间上保持零值不变，直到接下来的探测器单元发射的射线经过圆，才开始逐渐增加，当某一射线经过圆心时，再次达到一个峰值，随后减为零，边缘只穿过空气的射线衰减仍为零。

以图9为例，建立平面直角坐标系，其中直线ac为椭圆切线，直线od经过椭圆中心，直线be经过圆心，假设ac与od间的距离为x1，od与be间的距离为x2。设ac的直线方程为y＝k(x+a)，经过以上对接收信息的分析，可以发现x1与x2的比值为图中x3与x4的比值，再由简单的几何关系知

观察图9得椭圆中心与圆心的距离ob，而x3和x4的比值可由模板的接收信息确定，从而可求a0的长度，进而求出直线的横截距a1。

图5与图9是对称的情况

图6与图8也是对称的情况，参照分析图8时的几何方法，只是选取的相似三角形不一样，利用相似三角形对应边成比例的定理，得到图5、图7横截距满足的比例关系为

图7中由于射线方向与x轴方向平行，所以斜率为0。

把以上各情况下的横截距带入直线方程并与椭圆方程联立求解直线斜率。

y＝k(x+a1)

整理可得

(a²+a²k²)x²+2a²k²a1x+a²k²a1²-a²b²＝0

利用判别式法建立ct系统的旋转方向模型

δ＝2a²k²a1-4(a²+a²k²)(a²k²a1²-a²b²)＝0(4)

由图2的信息可求出椭圆方程的长轴长a＝15与短轴长b＝40，将a,b以及(3)式求解出的a1代入(4)式即可求出直线斜率，直线斜率即为直线与x轴正方向的夹角的正切值，利用matlab编程(见附录程序1)求得各直线斜率与旋转角度，求出旋转角度的初始值为60.3664度，求得的180个方向下的旋转角度见附录表1。

2)ct系统的探测器单元间距模型二

椭圆的另一条切线的确定需要其横截距的值，此时k是已知的，观察(4)式为关于横截距的一元二次方程，公式法即可求解。当两条平行直线的代数方程都确定时，两者间的距离可表示为

其中a2表示另一切线的横截距，由于射线是一组垂直于探测器平面的平行直线，因此求得的d是经过椭圆的射线路径对应的所有探测器单元的长度之和。d0是所求的探测器单元的间距，n为d对应的探测器单元数，得到ct系统的探测器单元间距模型二如下

利用matlab将数据分为如上五种类型，按照如上步骤的算法思想，将求得的结果取平均值得到间距为0.2753。

3)ct系统的旋转中心模型三

如图11所示，直线a和直线b分别为在探测器方向1、探测器方向2下，垂直于探测器平面的中垂线，直线c过原点，由于椭圆关于原点对称，所以直线c到直线a、b的距离相等。通过观察模板的接收信息的数据，可求出直线b、c之间对应的探测器单元总数，由于间距已经算出，因此b、c之间的距离以及整个探测器的长度可以求出。直线b、c之间的距离为

其中b0为直线b的总截距，由于各个方向的斜率已经求出，因此可以确定出直线b的斜截式，用同样的方法可以求出直线a的方程，联立两条直线方程

代入由(5)算出的横截距可求解出两条中垂线的交点

点(x0,y0)就是需要求得旋转中心的位置在图示坐标系下的位置。

为了避免偶然性，随机取模板的接收信息中的5列数据对应的射线方向，射线方向对应的斜率已求出，两两匹配求得的坐标取平均值得到旋转中心的坐标为(-9.4014,6.1567)。

二、未知介质1的位置、几何形状和吸收率的确定方法如下：

ct成像相当于对一切探测出的投影值rg(t,θ)，反求衰减系数的分布函数g(x,y)，这个反演过程可通过拉东变换实现；

拉东变换的参数形式定义式为

其中g(x,y)函数沿平面上任何一条直线xcosθ+ysinθ＝t可积。上式称为函数g(x,y)的拉东变换，记b(t,θ)＝rg(t,θ)；

先求出ct系统的增益系数，然后对数据进行去增益化

其中α＝ln[(l0-l)l]，计算出α＝1.75；

由于ct系统的旋转中心不在正方形托盘的中心，因此需要对数据做平移变化才能恢复数据的原始信息；

设旋转中心为(a,b)，做如下平移变化

rρ(t,θ)＝rf(t-acosθ-bsinθ,θ)

由中心切片定理知，图像重构就是求投影g(t,θ)的逆变换(反演变换)，对当θ给定时，图像投影rg(t,θ)的一维傅里叶变换等于对图像g(x,y)的二维傅里叶变换；

首先对图像投影做一维傅里叶变换

再做二维傅里叶变换

此时通过观测值b(t,θ)反推出的f(x,y)为经过校正后的信息。

通过matlab编程(见附录程序2)计算出图3中10个点的吸收率分别为:0.0001,0.0001.0.0003,0.0001,0.0002,0.0001,0.0001,0.0002,0.0004,1.000，呈现的图像如图12。

三、未知介质2的相关信息的确定方法如下：

未知介质2的位置、几何形状和吸收率的确定方法与未知介质1的相关信息的理论求解过程一致，相互区别的只是相应的经过增益后的数据不同。利用matalb自带的函数wdencmp函数进行小波去噪，重建出未知物质2的的图像如图13所示。求得的图3中的10个点的吸收率分别为0.0001，0.0001，1.0000，1.0000，0.5004，1.0000，1.0000，1.0000，0.0500，1.0000。

四、ct系统的误差分析与稳定性分析的具体方法如下：

在ct系统的旋转次数是有限的情况下，图像重建是不适定的，利用被测物体的先验信息对问题进行约束，可求出原问题的近似解，正则化算法就是这样一种方法；

正则化方法的求解过程为在对线性代数方程组ag＝b进行求解时，可将此问题转化为计算||ag-b||2的极小值，但对于不适定问题的求解，可添加一个惩罚项，此时原问题变为：

其中，ψ(g)取为正则函数，为附加在正则函数上的惩罚项因子，对极小值进行求解也就是对公式最优解进行求解；

正则化方法从本质来讲，属于最小二乘方法的一种：

由于上述ct系统的参数的给出是建立在把数据看成连续的图像，将点a视为一列数据的最大值，但是由于数据是离散的，该列数据中的最大值不一定就是a处的值，若以a所在的值对应数据中的最大值计算，则在问题一中用中垂线法求解，中垂线的方程携带了误差。

通过matlab的nlinfit函数对模板的接收信息进行非线性拟合，求出点a的横纵坐标，以a的横坐标代替原本通过数据的增减变化求得的横坐标可以减少该方法引起的误差。

中垂线与椭圆两条切线的平分线之间的距离为

由此求出直线的截距，而斜率已知，求得的旋转中心的坐标仍满足

利用中垂线法求出旋转中心的横坐标为7.3300，误差为0.224。

下面开始讨论稳定性，小圆半径对距离稳定性的影响。根据几何关系，可以计算出探测点元之间的距离

当斜率为定值时，x1/x2也为一定值，由此可以计算出r的系数为(x1+1/x2)，因为r和x1/x2代入算出距离相对较小，由此可以判断小圆半径对稳定性的影响较小。

符号说明：

附录

表1

程序1

a＝[304302301301298296295294292291291288287286285284282282281280279278278277277276276276276276276277277278278279280281282283284285287288290292293295297299301303305307309311313315317319322]；

b＝[219216216216216216217217217217217217217218218218218218219219219219220220220221221221222222222223223224224225225225227226227227228228229229230230231231232233233234234235236236237237238]；

c＝[77767473727170696767666563636362616160606060595959595959596060606161626263646465666768697071737475777879818384868890919395]；

x1＝abs(b-a)；

x2＝abs(c-b)；

x1/x2；

k＝-40./sqrt((45.*x1./x2).^2-225)；％x射线斜率

k＝atan(k)；％计算弧度

k.*180./pi％计算角度

m＝45.*x1./x2％y＝k(x+m)％计算图中的射线的λ值，即两条射线之间的距离(m+45).*k./sqrt(k.^2+1)

数据

a＝[304302301301298296295294292291291288287286285284282282281280279278278277277276276276276276276277277278278279280281282283284285287288290292293295297299301303305307309311313315317319322]；

b＝[219216216216216216217217217217217217217218218218218218219219219219220220220221221221222222222223223224224225225225227226227227228228229229230230231231232233233234234235236236237237238]；

c＝[77767473727170696767666563636362616160606060595959595959596060606161626263646465666768697071737475777879818384868890919395]；

x1./x2

vpa(k)

a-m

k＝40./sqrt((45.*x1./x2).^2-225)

a＝[188185184182179177175173171168166164162]；

b＝[274274273272272271270.5270269268268267266]；

c＝[416413412409407405403400398395393391388]；

n-bh

k＝40./sqrt((45+45.*(x1-x2)./x2).^2-225)

a＝[160157155153151149147145143141139137135133131129127126124122120119117116114113111110108107105104103102100999897969594939392919090]；

b＝[266265264264263262262261260259.5259258257257256255255254253252252251250250249248248247246246246244243.5243242241.5241240240239238238238234231227224]；

c＝[385383380377363361359357354352350347345342339337334331329325324319316313310307304301298295291288285282278275272268265262258255252248245242238]；bi

bj-df

k＝-40./sqrt((45+45.*(x1-x2)./x2).^2-225)

a＝[378378377377376375374374373372371370369368367366365364363362361359358357356354353352351349348346345344342341339338336335333332330329327325324322321]；

b＝[246243240236233231230230229229228228227227226226225225224224223223223222222222221221221220220220219219219219218218218218218217217217217217217216216]；

c＝[2322282252222192152122092062021991961931901871841811781751721691661631601581551521501471451421401371351331311291271251231211191181161151131129694]；

dg-fx

k＝-40./sqrt((45.*x1./x2).^2-225)

a＝[319318316314313311310308306305303302300299296295294292291291288287286285284282282281280279278278277277276276276276276276277277278278279280281282283284285287288290292293295297299301303305307309311313315317319322]；

b＝[216216216216216216216216216216216216216216216217217217217217217217218218218218218219219219219220220220221221221222222222223223224224225225225227226227227228228229229230230231231232233233234234235236236237237238]；

c＝[92918987858382807977767473727170696767666563636362616160606060595959595959596060606161626263646465666768697071737475777879818384868890919395]；

程序2

d＝0.2768；

xc＝-33.5*d；

yc＝20*d；％旋转中心在托盘上的位置

phantom＝a；

phantom＝[zeros(100,180)；phantom；zeros(100,180)]；

img＝iradon(phantom,[0:179]+30)％通过拉冬变换求得吸收率

n＝size(img,1)；

[x,y]＝meshgrid([-n/2:n/2]*d)；％100/256imagesc(x(1,:),y(:,1),img)％重建影象

colormap(gray)％灰度

holdon

line([-5050],[-50-50],'color','w')

line([-5050],[5050],'color','w')

line([5050],[-5050],'color','w')

line([-50-50],[-5050],'color','w')％托盘位置

holdon

plot(xc,yc,'ok')

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。