专利名称：工业ct系统探测器分组串绕的校正方法
技术领域：
本发明属在核技术领域中用数学方法对CT系统重建图象进行校正。
2)对于解卷积迭代方法，因为要进行迭代，而且一般要经过多次迭代才能达到校正的目的，对于大数据量的工业CT系统并不实用。而且因为解卷积过程本身是病态过程，不能保证收敛，所以校正结果没有保证。
3)在杂志IEEE Trans.Image Processing，1995，4(10)1430-1438和IEEE Trans.Med.Imaging，1998，17(6)1101-1105.中分别刊载的题目为“Globally Convergent Algorithms for Maximum a PosterioriTransmission Tomography”，“Accelerated Iterative TransmissionCT Reconstruction Using an Ordered Subsets Convex Algorithm”的文献分别提供了“Convex方法”及“Ordered Subset Convex方法”(序列子集凸函数方法)及其标准过程，Convex方法以从源发出的光子数和探测器接收到的光子数符合泊松分布为基础，得出探测器接收光子数的对数似然函数，当此对数似然函数取得极大值时所得的输入吸收系数uj值即是重建所述的uj值，uj通过迭代过程求出；OrderedSubsets Convex方法是将投影数据分为若干子集，对每个子集执行一次Convex迭代，所得结果作为下一个子集迭代的输入，处理完所有子集就是完成了一次Ordered Subsets Convex迭代。
4)对于迭代重建中加入探测器串绕模型的方法，目前只限于对一般探测器串绕的校正，而且校正速度慢，需要多次迭代才能完成。而且目前的迭代方法中很多方法收敛很慢，对于实际的系统不适用。
本发明的技术方案如下以序列子集凸函数方法即Ordered Subsets Convex方法为基础，用探测器分组串绕模型处理每次迭代产生的投影数据，并在每次迭代后进行中值滤波，具体步骤如下(1)以未校正的重建图象作为迭代开始的初始吸收系数uj值输入；
(2)将实际得到的探测器数据分为几个子集，从第一个子集开始处理；(3)用输入图象计算积分值<lit，u>；(4)将积分值求指数后进行抽取，卷积，再求对数后放回原位置；(5)计算新的输入吸收系数u值，求得新的重建图象；(6)对所得新重建图象进行中值滤波；(7)所有子集处理完了没有？若还没处理完，则以滤波后的图象作为新的输入图象返回步骤(3)，对下一个子集进计算，若已全部处理完，则进入步骤(8)；(8)检验迭代次数是否达到预定迭代次数？如果没有达到预定迭代次数，则以滤波后的图象作为输入图象返回步骤(2)继续迭代；(9)已达到预定迭代次数，停止处理。
本发明的积极效果是1、由于在每次迭代后对得到新的重建图象进行中值滤波加快了收敛，因此所需迭代次数少，一般情况一次迭代即可完成，对于串绕严重的情况，三次迭代以内也可完成校正，并不显著增加计算时间；2、由于采用了Ordered Subsets Convex迭代，并用探测器分组串绕模型加入迭代重建过程，所以校正质量好，校正完成后所得图象的灰度值和待重建物体的吸收系数对应好，达到了很好的校正目的。
图中·为校正后数据，*为校正前数据，+为真实数据，比较所用数据为待重建截面中横坐标为46的截线段。
对数技术方案中的步骤(5)所述的求得新的重建图象，具体方法是用步骤(4)中所得新的积分值带入ujn+1=ujnΣiΣtlitj[dite-<lit,un>(1+<lit,un>)-Yit]ΣiΣtlitj<lit,un>dite-<lit,un>]]>计算得出新的重建图象在每一点的灰度值，其中i＝(s-1)*r/n+1，…，s*r/n，t＝1，2，…，k，Yit是在i投影方向在探测器t处实际接收到的光子数，dit是在投影方向i时沿着探测器t方向从源发出的光子数。
以下结合图3、图4对本发明技术方案的实施作进一步说明。
设扫描系统参数为每次扫描由k＝128个探测器数据组成，探测器每8个形成一个探测器组，同组探测器发生串绕，源绕物体一周共进行r＝360次扫描。待重建图象有m＝128*128＝16384个象素点。子集数取n＝36。探测器串绕函数取5点离散高斯函数 ，中值滤波模板取图3所示二维5点滤波模板，迭代2次后结束。具体步骤为1)以未校正的重建图象作为迭代开始的初始吸收系数uj(j＝1，2，Λ，16384)。
2)将实际得到的探测器数据按照扫描时从0度到360度的旋转平均分为36个子集，则每个子集包含10个投影方向。从第一个子集开始按照下述步骤计算3)设子集数为s，在某一投影方向i上，求象素j在t探测器方向上的相交长度litj。然后得到输入吸收系数uj在各个探测器方向t上的积分值<lit,u>=Σjlitjuj]]>。对于t＝1，2，…，128，计算出所有积分值<lit，u>。对于该子集内的所有i，即i＝(s-1)*10+1，…，s*10，也重复上述步骤，计算出所有的积分值<lit，u>(i＝(s-1)*10+1，…，s*10，t＝1，2，…，128)。然后对所有积分值<lit，u>取指数得到投影值pit=e-<lit,u>.]]>4)将相同投影方向i时的投影值pit按下标t顺序排列，得到投影值列向量Pi=pi1pi2Mpi128]]>(i＝(s-1)*10+1，…，s*10)，对得到的投影值矩阵Pi中的投影值Pit，抽出Pi1，Pi9，Pi17…Pi121为一组，Pi2，Pi10，Pi18…Pi122为一组，Pi3，Pi11，Pi19…Pi123为一组，依次类推，每个Pi得到16组。每组投影值都和探测器串绕函数f(x)进行卷积，用所得卷积结果代替原来的投影值Pit，并用-ln(Pit)作为新的积分值<lit，u>。各组数据都处理完后仍放回抽取前的位置。
5)用步骤(4)中所得新的积分值带入ujn+1=ujnΣiΣtlitj[dite-<lit,un>(1+<lit,un>)-Yit]ΣiΣtlitj<lit,un>dite-<lit,un>,(i=(s-1)*10+1,···,]]>s*10，t＝1，2，…，128)，计算得出新的重建图象在每一点的灰度值。
6)对所得新重建图象用图3所示模板进行中值滤波。
7)返回步骤(3)对下一个子集s＝s+1进行计算，直到处理完所有子集s＝36。
8)返回步骤(2)继续迭代，如果已经是第2次迭代则停止计算。
权利要求
1.一种工业CT系统探测器分组串绕的校正方法，包括序列子集凸函数方法即Ordered Subsets Convex方法，探测器串绕模型加入迭代重建，以及二维中值滤波，其特征在于以序列子集凸函数方法即Ordered SubsetsConvex方法为基础，用探测器分组串绕模型处理每次迭代产生的投影数据，并在每次迭代后进行中值滤波；具体步骤如下(1)以未校正的重建图象作为迭代开始的初始吸收系数值输入；(2)将实际得到的探测器数据分为几个子集，从第一个子集开始处理；(3)用输入图象计算积分值<lit，u>；(4)将积分值求指数后进行抽取，卷积，再求对数后放回原位置；(5)计算新的输入吸收系数u值，求得新的重建图像；(6)对所得新重建图象进行中值滤波；(7)所有子集处理完了没有？若还没处理完，则以滤波后的图象作为新的输入图象返回步骤(3)，对下一个子集进行计算，若已全部处理完，则进入步骤(8)；(8)检验迭代次数是否达到预定迭代次数，如果没有达到预定迭代次数，则以滤波后的图象作为输入图象返回步骤(2)继续迭代；(9)已达到预定迭代次数，停止处理。
2.根据权利要求1所述的工业CT系统探测器分组串绕的校正方法，其特征在于步骤(3)所述用输入图象计算积分值的具体方法是设子集数为S，在某一个投影方向上，以输入吸收系数u计算在各个探测器方向t上的积分值<lit,u>=Σjlitjuj]]>，其中litj是象素j在i投影方向时在t探测器方向上的相交长度；对t计算出所有的积分值，<lit，u>，对于该子集内的所有i重复上述步骤，计算出所有的积分值<lit，u>。
3.根据权利要求1或2所述的工业CT系统探测器分组串绕的校正方法，其特征在于所述相交长度litj可以采用近似计算得到，具体方法是如果投影方向i时从源到探测器t的射线不经过点j，则litj＝0，如果经过，设相交长度litj与X轴所夹的锐角为θ，以X轴和Y轴的单位长度作为1，当相交长度litj与X轴的夹角θ小于45度时litj＝1/cos(θ)，当夹角大于45度时litj＝1/sin(θ)；点j处的输入吸收系数uj与对应的相交长度litj相乘，并对下标j求和，得<lit,u>=Σjlitjuj.]]>
4.根据权利要求1所述的工业CT系统探测器分组串绕的校正方法，其特征在于步骤(4)所述将积分值求指数后进行抽取，卷积，再求对数后放回原位置的具体方法是对步骤(4)中得到的所有积分值<lit，u>取指数得到投影值pit=e-<lit,u>]]>；将相同投影方向i时的投影值pit按下标t的顺序排列，则每个投影方向上都得到投影值列向量Pi=pilpi2Mpiki]]>对得到的投影值矩阵Pi中的投影值pit，按照实际的探测器分组间隔进行抽取，得到不同的投影值组，每组投影值和探测器串绕函数f(x)进行卷积，用所得卷积结果代替原来的投影值Pit，并取该投影值的自然对数的负值-ln(pit)作为新的积分值<lit，u>，各组数据都处理完后仍放回抽取前的位置。
5.根据权利要求4所述的工业CT的系统探测器分组串绕的校正方法，其特征在于所述探测器串绕函数f(x)可以从实验中测得，将射线束对着一组探测器的中心探测器，x坐标轴中心置于该中心探测器，x方向沿探测器分布方向，以射入射线的强度作为1，测量该中心探测器，以及其它相临探测器的读数，所得不同x位置处的读数即是探测器串绕函数f(x)。
6.根据权利要求1所述的工业CT系统探测器分组串绕的校正方法，其特征在于步骤(5)所述的求得新的重建图像，具体方法是用步骤(4)中所得新的积分值带入ujn+1=ujnΣiΣtlitj[dite-<lit,un>(1+<lit,un>)-Yit]ΣiΣtlitj<lit,un>dite-<lit,un>]]>计算得出新的重建图象在每一点的灰度值，其中Yit，是在i投影方向在探测器t处实际接收到的光子数，dit是在投影方向i时沿着探测器t方向从源发出的光子数。
全文摘要
工业CT系统探测器分组串绕的校正方法，属在核技术领域中用数学方法对CT系统重建图象进行校正，主要特征是以序列子集凸函数方法Ordered Subsets Convex为基础，用探测器分组串绕模型处理每次迭代产生的投影数据，并在每次迭代后进行中值滤波以加快收敛。本发明的校正方法所需迭代次数少，校正质量好，适于实际的工业CT系统中重建图象的校正，克服了传统技术中迭代次数多，校正效果不好等问题。
文档编号G01N23/02GK1407439SQ0113069
公开日2003年4月2日 申请日期2001年8月30日 优先权日2001年8月30日
发明者孙少华, 高文焕, 张丽, 陈志强 申请人:清华大学, 清华同方威视技术股份有限公司