CN108364325A - 规则样品x射线ct投影图位置平移偏差检测及校正方法 - Google Patents

Abstract

本发明涉及CT成像技术领域，特别涉及规则样品X光断层扫描投影图像中样品发生随机平移后对偏移量的检测及其校正方法。在原始CT投影图或者减背景投影图的样品区域中选取判断样品偏移量的参考点；计算每张投影图对应位置处样品发生的偏移量，并进行平移修正。本方法结合了三个不定参数进行拟合，并完成对样品发生随机平移的图像的校正，以此减少样品随机偏移造成的重构数据中组分的形状畸变，为后续的CT重构提供良好的数据基础。

Description

规则样品X射线CT投影图位置平移偏差检测及校正方法

技术领域

本发明涉及CT成像技术领域，特别涉及规则样品X光断层扫描投影图像中样品发生随机平移后对偏移量的检测及其校正方法。

背景技术

计算机断层扫描(CT，Computed Tomography)技术可以不破坏样品结构而获取样品内部结构信息，具有无损性、高分辨率、高灵敏度等成像技术特点。被广泛应用于医学、工业、生物、材料科学等领域。近年，同步辐射X射线CT技术较普通X射线管光源具有优良的单色性，被广泛应用于前沿科学研究，成为一种重要的材料结构表征技术。

CT成像质量是影响材料结构表征精确性的重要因素，高精度的CT图像重建要求样品在旋转扫描过程中，光源、转台中心及探测器保持相对静止，否则在重建的图像中将出现组分形状畸变，严重影响重构图像的质量。但是由于机械振动、机械精度以及环境因素的干扰，在扫描过程中光源、转台或者背景总是会发生随机偏移。

目前，有很多关于CT转台的旋转中心随机偏移的校正方法，常用的方法有采用特殊设计的校正模体进行校正；针对CT投影图中样品位置发生偏移，常用的方法有利用第一张投影图和最后一张投影图的偏差量，然后将总的偏差量平均分配到每一张投影图中，再对每一张图片进行相应的平移，但是此方法适用于样品旋转360°的情况。这些方法在特定的CT扫描系统上都获得了较好的效果。但是，在CT扫描过程中，不仅旋转中心随机偏移，也有可能是样品本身发生偏移。本发明提出一种在投影图中基于样品参考点坐标的多参数最优化校正技术，不管是样品偏移还是样品台偏移引起投影图中样品位置发生偏移，此方法都可以对其进行修正，只需由采集的投影图片来进行修正，适用于样品台旋转任何角度的情况，且实现起来简单方便，实用性较强。

发明内容

本发明所要解决的技术问题是：如何在规则样品X光断层扫描投影图像中样品发生随机平移后对偏移量进行检测及其校正，以获得较好的重建图像。

本发明所采用的技术方案是：规则样品X射线CT投影图位置平移偏差检测及校正方法，按照如下的步骤进行

步骤一、以采集的CT投影图上左上角的第一个像素点为坐标原点建立直角平面坐标系，向右表示投影图横坐标i的正方向，i＝0,1,2,…i_max，i_max是投影图的最大横坐标值；向下表示图像纵坐标j的正方向，j＝0,1,2,…j_max，j_max是投影图的最大纵坐标值， (i_max+1)×(j_max+1)为投影图的像素尺寸，在CT投影图的样品区域中选取一个有标志性的点作为判断样品偏移量的参考点，确定其实验横坐标值直角平面坐标系，以一个像素点为最小单位，本专利所述的规则样品是指样品的几何图形为规则状的长方体、圆柱体、棱锥体、球体等符合几何规则的图形，不规则的图形不适合本专利，注意，本专利中针对的是任意一张采集的CT投影图，其它张的投影图的处理方法都与本专利类似；

步骤二、在步骤一所建立的直角平面坐标系中，建立CT投影图中参考点的理论横坐标值的函数这一函数中包括样品台旋转中心在投影图中的投影横坐标i_s、参考点离样品台旋转中心的距离r、规则样品转动初始相位θ三个参数；其中各自的取值范围为0＜i_s＜i_max、0≤r＜i_max/2、0≤θ≤2π；k表示第k张CT投影图，φ为样品台转动过程中的角度步长；

步骤三、建立确定参考点横坐标理论值函数所包含变量最优解的目标函数，

k_max表示CT投影图的总张数；

步骤四、通过使得目标函数取到极小值，以此确定参考点横坐标的理论值函数中三个参数： i_s、r、θ的最优解；

步骤五、将i_s、r、θ的最优解代入步骤二所建立的参考点理论横坐标值的函数中，计算每张投影图中参考点的理论横坐标值再根据参考点的实验横坐标值计算出每张CT 投影图偏移量表示取数值的整数部分；

步骤六、根据偏移量x^k，对第k张CT投影图进行平移校正，本专利中认为偏差的产生主要是由于样品旋转台、规则样品在转动过程中发生震动或者位移而产生的，因此每张CT投影图的每个像素点偏差都认为完全相同；

步骤七、用同样的方法对所有CT投影图像位置平移偏差进行检测及校正。

作为一种优选方式：所述步骤一，在CT投影图的样品区域中选取一个有标志性的点作为判断样品偏移量的参考点，确定其实验横坐标值其具体过程为，读取原始CT投影图像的灰度值，建立灰度数值矩阵I_k(i,j)，表示第k张投影图像中坐标为(i,j)的像素点的灰度值；在CT投影图中所建立的灰度数据矩阵中样品所在区域内选取n行1≤n≤j_max，计算其每一列像素点的平均灰度值，然后计算由n行得到的每一列平均灰度值的样品左右边界实验横坐标，其实验纵坐标是所选取n行中的任意一行，设纵坐标为j_n′，n行中第一行的纵坐标值为j_m，n行中最后一行的纵坐标值为j_m+n-1，0≤j_m＜j_m+n-1≤j_max；对第j_n′行的区间 [Δi,i_max-Δi]中每个点的平均灰度值进行差分计算，得到左右差分函数的乘积f_k′(i)，并由其最小值计算样品边界点的实验横坐标：

点i的取值范围为(Δi≤i≤(i_max-Δi))，其中Δi为大于0的正整数，要保证该行样品边界点的 f_k′(i)大于噪声涨落最大点的f_k′(i)的前提下取最小正整数，噪声涨落最大点是由投影图背景噪声情况来确定，I_k(i,j)表示第k张投影图像中实验坐标为(i,j)像素点的灰度值；在区间Δi≤i≤(i_max+1)/2内，当f_k′(i)取到最小值f_k′(i₁)_min时，对应自变量为样品左边界坐标；在区间(i_max+1)/2＜i≤(i_max-Δi)内，当f_k′(i)取到最小值f_k′(i₂)_min时，对应自变量为样品右边界坐标；即与为第k张CT投影图第j_n′行平均灰度值的样品左右边界实验横坐标，则该行样品中心横坐标的实验值为并以此作为判断偏移量的参考点。

作为一种优选方式：所述步骤四中，通过使得目标函数取到极小值，以此确定参考点横坐标的理论值函数中三个参数：i_s、r、θ的最优解；

方法1、首先选取i_s＝i_max/2、r＝0、θ＝0的初始值作为起点，然后根据数学中的常规方法爬山法确定目标函数的极小值，从而确定三个参数i_s、r、θ的最优解，在利用爬山法计算目标函数极小值的过程中，三个参数i_s、r、θ的步长在范围0＜Δi_s＜0.5、0＜Δr＜1 与0＜Δθ≤φ中可以取任意确定值。

方法2、首先选取i_s＝i_max/2、r＝0、θ＝0的初始值作为起点，然后找到Y(i_s+Δi_s,r,θ)、 Y(i_s-Δi_s,r,θ)、Y(i_s,r+Δr,θ)、Y(i_s,r-Δr,θ)、Y(i_s,r,θ+Δθ)、与Y(i_s,r,θ-Δθ)六个函数值中的极小值并相应的改变其中的参数，然后将这个极小值依据上面的方法再构造一个解集，再次求得极小值，以此类推，直到后一次的极小值比前一次的极小值大就结束，其中Δi_s、Δr、Δθ分别表示三个参数的步长，取值范围分别为0＜Δi_s＜0.5、0＜Δr＜1与0＜Δθ≤φ，当三个参数的步长在规定的范围取一定值，由目标函数的极小值计算得到i_s、r、θ三个参数的近似最优解作为最优解。

需要注意的是，本发明提供了两种获得最优解的方法，但是不是说就是只有这两种方法，实际上求最优解的过程是一个数学方法，可能还存在其它的求得最优解的方法。

作为一种优选方式：所述步骤六中，根据偏移量x^k，对第k张CT投影图进行平移校正的具体过程为，I_k(i,j)表示第k张投影图像中实验坐标为(i,j)像素点的灰度值，平移修正后，I_k′(i′,j′)表示第k张修正后投影图像中理论坐标为(i′,j′)像素点的灰度值，且i′的取值范围为0≤i′≤i_max，j′的取值范围为0≤j′≤j_max，

当x^k≥0时：

当x^k＜0时：

本发明的有益效果是：通过本发明方法使修正后的投影图重构的CT切片中的组分的形状畸变得到修复，样品细节得到准确复现。

附图说明

图1为页岩的CT投影图；

图2为页岩减背景后的投影图；

图3为分析参数i_s,r,θ三者关系的简易图；

图4为截取平移修正前的CT投影图重建的CT切片图；

图5为所截取的平移修正后的CT投影图重建的CT切片图。

具体实施方式

本发明所述的规则样品X射线CT投影图位置平移偏差检测及校正方法，可用来修正 CT设备获取的规则样品投影图像，下面将结合附图对本发明进行详细描述。

本实施例所采用投影数据是在上海光源BL13W线站采集的页岩CT投影图，样品直径为1.0mm，分辨率为0.65um，投影图像素尺寸为2048×2048，本次实验样品台旋转180°共采集1260张投影图。本实施例只是为了进一步阐述本发明，但并不限制本发明所保护的范围。

规则样品X射线CT投影图位置平移偏差检测及校正方法，包括以下步骤：

步骤一、以采集的CT投影图上左上角的第一个像素点为坐标原点建立直角平面坐标系，向右表示投影图横坐标i的正方向，i＝0,1,2,…i_max，i_max是投影图的最大横坐标值，本实施例中i_max＝2047；向下表示图像纵坐标j的正方向，j＝0,1,2,…j_max，j_max是投影图的最大纵坐标值，本实施例中j_max＝2047，由于采集的CT投影图背景噪声影响比较大，所以首先进行了减背景处理如图2，然后读取减背景CT投影图像的灰度值，建立灰度数值矩阵I_k(i,j)，表示第k张投影图像中坐标为(i,j)的像素点的灰度值；在减背景CT投影图中所建立的灰度数据矩阵中样品所在区域内选取n行1≤n≤j_max，计算其每一列像素点的平均灰度值，然后计算由n行得到的每一列平均灰度值的样品左右边界实验横坐标，其实验纵坐标是所选取n 行中的任意一行，设纵坐标为j_n′，n行中第一行的纵坐标值为j_m，n行中最后一行的纵坐标值为j_m+n-1，0≤j_m＜j_m+n-1≤j_max；对第j_n′行的区间[Δi,i_max-Δi]中每个点的平均灰度值进行差分计算，得到左右差分函数的乘积f_k′(i)，并由其最小值计算样品边界点的实验横坐标：

点i的取值范围为(Δi≤i≤(i_max-Δi))，其中Δi为大于0的正整数，要保证该行样品边界点的 f_k′(i)大于噪声涨落最大点的f_k′(i)的前提下取最小正整数，噪声涨落最大点是由投影图背景噪声情况来确定，I_k(i,j)表示第k张投影图像中实验坐标为(i,j)像素点的灰度值；在区间Δi≤i≤(i_max+1)/2内，当f_k′(i)取到最小值f_k′(i₁)_min时，对应自变量为样品左边界坐标；在区间(i_max+1)/2＜i≤(i_max-Δi)内，当f_k′(i)取到最小值f_k′(i₂)_min时，对应自变量为样品右边界坐标；即与为第k张CT投影图第j_n′行平均灰度值的样品左右边界实验横坐标，则该行样品中心横坐标的实验值为并以此作为判断偏移量的参考点。

这里以选择第500张投影图进行说明；选取第500张投影图的900-915行(即n＝16)对每一列进行平均计算，每行宽度为1像素。当Δi＝1时，差分涨落参差不齐，很难找到f₅′₀₀(i) 的极值；当Δi＝2时，除样品边界点外其它点对应的f₅′₀₀(i)都涨落比较均匀，选择Δi＝2较为理想。由公式计算得到在2≤i≤1024的区间内，差分函数最小值f₅′₀₀(i₁)_min＝-15192.3，相应的为左边界；在1024＜i≤2045的区间内，差分函数最小值 f₅′₀₀(i₂)_min＝-5324.37，相应的为右边界；利用公式计算得到样品的中心横坐标

步骤二、在步骤一所建立的直角平面坐标系中，依据对所有CT投影图中参考点横坐标的实验值进行拟合，建立CT投影图中参考点的理论横坐标值的函数这一函数中包括样品台旋转中心在投影图中的投影横坐标i_s、参考点离样品台旋转中心的距离r、样品转动初始相位θ三个参数；其中各自的取值范围为0＜i_s＜i_max、0≤r＜i_max/2、 0≤θ≤2π；k表示第k张CT投影图，φ为样品台转动过程中的角度步长；如图3所示，P 点为样品，O点为投影图的零坐标，其中表示第k张投影图中参考点横坐标的理论值 |OO₃|，i_s为样品台旋转中心在CT投影图的投影横坐标|OO₂|；r为参考点离样品台旋转中心的距离|PO₁|；φ为样品台转动过程中的角度步长，实施例中φ＝π/1260；θ为样品转动的初始相位，范围为(0≤θ≤π)。

步骤三、建立确定参考点横坐标理论值函数所包含变量最优解的目标函数：

Y(i_s,r,θ)是为确定样品参考点横坐标理论值函数其所包含的三个变量i_s、r、θ的最优解所建立的目标函数，本实施例中k_max＝1260，本专利中Y(is,r,θ)简写为Y。

步骤四、通过使得目标函数取到极小值，以此确定参考点横坐标的理论值函数中三个参数： i_s、r、θ的最优解；

方法1、首先选取i_s＝i_max/2、r＝0、θ＝0的初始值作为起点，然后根据数学中的常规方法爬山法确定目标函数的极小值，从而确定三个参数i_s、r、θ的最优解，在利用爬山法计算目标函数极小值的过程中，三个参数i_s、r、θ的步长在范围0＜Δi_s＜0.5、0＜Δr＜1与0＜Δθ≤φ中可以取任意确定值。

方法2、首先选取i_s＝i_max/2、r＝0、θ＝0的初始值作为起点，然后找到Y(i_s+Δi_s,r,θ)、 Y(i_s-Δi_s,r,θ)、Y(i_s,r+Δr,θ)、Y(i_s,r-Δr,θ)、Y(i_s,r,θ+Δθ)、与Y(i_s,r,θ-Δθ)六个函数值中的极小值并相应的改变其中的参数，然后将这个极小值依据上面的方法再构造一个解集，再次求得极小值，以此类推，直到后一次的极小值比前一次的极小值大就结束，其中Δi_s、Δr、Δθ分别表示三个参数的步长，取值范围分别为0＜Δi_s＜0.5、0＜Δr＜1与0＜Δθ≤φ，当三个参数的步长在规定的范围取一定值，由目标函数的极小值计算得到i_s、r、θ三个参数的近似最优解作为最优解。

需要注意的是，本发明提供了两种获得最优解的方法，但是不是说就是只有这两种方法，实际上求最优解的过程是一个数学方法，可能还存在其它的求得最优解的方法。

计算最优解时选取三个参数的步长分别为Δi_s＝0.1、Δr＝0.1与Δθ＝φ，计算得到 Y_min(s,r,θ)＝22788，对应得到三个参数的最优解s＝1021.8,r＝118.8(单位为像素)；θ＝65.8379°。

步骤五、将i_s、r、θ的最优解代入步骤二所建立的参考点理论横坐标值的函数中，计算每张投影图中参考点的理论横坐标值再根据参考点的实验横坐标值计算出每张CT 投影图偏移量表示取数值的整数部分；

当k＝500时，由步骤二得到第500张图片样品参考点坐标的理论值然后再依据其实验横坐标计算得到第500张图片样品的偏移量：

x⁵⁰⁰＝floor(1090-1098.862+0.5)＝-9，单位为像素。

步骤六、根据偏移量x^k，对第k张CT投影图进行平移校正的具体过程为，I_k(i,j)表示第k 张投影图像中实验坐标为(i,j)像素点的灰度值，平移修正后，I_k′(i′,j′)表示第k张修正后投影图像中理论坐标为(i′,j′)像素点的灰度值，且i′的取值范围为0≤i′≤i_max，j′的取值范围为0≤j′≤j_max，

当x^k≥0时：

当x^k＜0时：

由上计算的第500张CT投影图对应的偏移量x⁵⁰⁰＝-9，单位为像素；由于x⁵⁰⁰＜0，由以上公式得到第500张减背景后的投影图向右平移9个像素。

步骤七、利用计算机以同样方法对所有CT投影图像的平移偏差进行检测及其校正，并输出这些图像。图4为平移修正前的CT投影图重建的CT切片图，图5为平移修正后的CT投影图重建的CT切片图。通过对比图4与图5可以看出，修正后的投影图重构的CT切片中的组分的形状畸变明显减少了，样品细节得到准确复现。

Claims (4)

1.规则样品X射线CT投影图位置平移偏差检测及校正方法，其特征在于：按照如下的步骤进行

步骤一、以采集的CT投影图上左上角的第一个像素点为坐标原点建立直角平面坐标系，向右表示投影图横坐标i的正方向，i＝0,1,2,…i_max，i_max是投影图的最大横坐标值；向下表示图像纵坐标j的正方向，j＝0,1,2,…j_max，j_max是投影图的最大纵坐标值，

(i_max+1)×(j_max+1)为投影图的像素尺寸，在CT投影图的规则样品区域中选取一个有标志性的点作为判断规则样品规则偏移量的参考点，确定其实验横坐标值

步骤二、在步骤一所建立的直角平面坐标系中，建立CT投影图中参考点的理论横坐标值的函数这一函数中包括样品台旋转中心在投影图中的投影横坐标i_s、参考点离样品台旋转中心的距离r、规则样品转动初始相位θ三个参数；其中各自的取值范围为0＜i_s＜i_max、0≤r＜i_max/2、0≤θ≤2π；k表示第k张CT投影图，φ为样品台转动过程中的角度步长；

步骤三、建立确定参考点横坐标理论值函数所包含变量最优解的目标函数，

k_max表示CT投影图的总张数；

步骤四、通过使得目标函数取到极小值，以此确定参考点横坐标的理论值函数中三个参数：i_s、r、θ的最优解；

步骤五、将i_s、r、θ的最优解代入步骤二所建立的参考点理论横坐标值的函数中，计算每张投影图中参考点的理论横坐标值再根据参考点的实验横坐标值计算出每张CT投影图偏移量表示取数值的整数部分；

步骤六、根据偏移量x^k，对第k张CT投影图进行平移校正；

步骤七、用同样的方法对所有CT投影图像位置平移偏差进行检测及校正。

2.根据权利要求1所述的规则样品X射线CT投影图位置平移偏差检测及校正方法，其特征在于：所述步骤一，在CT投影图的规则样品区域中选取一个有标志性的点作为判断规则样品偏移量的参考点，确定其实验横坐标值其具体过程为，读取原始CT投影图像的灰度值，建立灰度数值矩阵I_k(i,j)，在CT投影图中所建立的灰度数据矩阵中规则样品所在区域内选取n行1≤n≤j_max，计算其每一列像素点的平均灰度值，然后计算由n行得到的每一列平均灰度值的规则样品左右边界实验横坐标，其实验纵坐标是所选取n行中的任意一行，设纵坐标为j_n′，n行中第一行的纵坐标值为j_m，n行中最后一行的纵坐标值为j_m+n-1，0≤j_m＜j_m+n-1≤j_max；对第j_n′行的区间[Δi,i_max-Δi]中每个点的平均灰度值进行差分计算，得到左右差分函数的乘积f_k′(i)，并由其最小值计算规则样品边界点的实验横坐标：

点i的取值范围为(Δi≤i≤(i_max-Δi))，其中Δi为大于0的正整数，要保证该行样品边界点的f_k′(i)大于噪声涨落最大点的f_k′(i)的前提下取最小正整数，噪声涨落最大点是由投影图背景噪声情况来确定，I_k(i,j)表示第k张投影图像中实验坐标为(i,j)像素点的灰度值；在区间Δi≤i≤(i_max+1)/2内，当f_k′(i)取到最小值f_k′(i₁)_min时，对应自变量为规则样品左边界坐标；在区间(i_max+1)/2＜i≤(i_max-Δi)内，当f_k′(i)取到最小值f_k′(i₂)_min时，对应自变量为规则样品右边界坐标；即与为第k张CT投影图第j_n′行平均灰度值的规则样品左右边界实验横坐标，则该行规则样品中心横坐标的实验值为并以此作为判断偏移量的参考点。

3.根据权利要求1所述的规则样品X射线CT投影图位置平移偏差检测及校正方法，其特征在于：所述步骤四中，通过使得目标函数取到极小值，以此确定参考点横坐标的理论值函数中三个参数：i_s、r、θ的最优解；首先选取i_s＝i_max/2、r＝0、θ＝0的初始值作为起点，然后找到Y(i_s+Δi_s,r,θ)、Y(i_s-Δi_s,r,θ)、Y(i_s,r+Δr,θ)、Y(i_s,r-Δr,θ)、Y(i_s,r,θ+Δθ)、与Y(i_s,r,θ-Δθ)六个函数值中的极小值并相应的改变其中的参数，然后将这个极小值依据上面的方法再构造一个解集，再次求得极小值，以此类推，直到后一次的极小值比前一次的极小值大结束，其中Δi_s、Δr、Δθ分别表示三个参数的步长，取值范围分别为0＜Δi_s＜0.5、0＜Δr＜1与0＜Δθ≤φ，当三个参数的步长在规定的范围取一定值，由目标函数的极小值计算得到i_s、r、θ三个参数的近似最优解作为最优解。

4.根据权利要求1所述的规则样品X射线CT投影图位置平移偏差检测及校正方法，其特征在于：所述步骤六中，根据偏移量x^k，对第k张CT投影图进行平移校正的具体过程为，I_k(i,j)表示第k张投影图像中实验坐标为(i,j)像素点的灰度值，平移修正后，I_k′(i′,j′)表示第k张修正后投影图像中理论坐标为(i′,j′)像素点的灰度值，且i′的取值范围为0≤i′≤i_max，j′的取值范围为0≤j′≤j_max，

当x^k≥0时：

当x^k＜0时：