CN103218834B - 一种基于点扩展函数的工业ct图像重建中心定位方法 - Google Patents

Abstract

一种基于点扩展函数的工业CT图像重建中心定位方法，涉及CT扫描图像构建领域，具体步骤如下：1）获取待测圆盘厚度方向中间部位断层投影数据，得到投影正弦图；2）计算边缘的质心位置x₀；3）设圆盘旋转中心偏移量x为x₀的±δ领域，依次运用滤波反投影重建算法重建各个偏移量x对应的CT图像；4）求出重建后每一幅CT图像中圆盘边缘的点扩展函数PSF的高度h，拟合图像重建中心偏移量x与中圆盘边缘的点扩展函数PSF的高度h之间的函数关系，计算图像重建中心偏移量。

Description

一种基于点扩展函数的工业CT图像重建中心定位方法

技术领域

本发明涉及CT扫描图像构建领域，特别是一种用于工业CT图像重建中心定位的方法。

背景技术

Ⅲ代CT扫描图像重建是基于这样一种假设，扇形射束在旋转投影时，中心射线经过旋转中心，或者说旋转中心在扇形射线束中的位置是已知的，实际CT系统要做到这一点是很困难的，在设备最初调试阶段要经过繁琐的试凑过程才能满足精度需求；而且当设备使用一段时间，因环境和负载变化等因素，旋转中心位置会发生漂移，找准位置费力费时。针对这一问题解决方法很多，主要有迭代法，正弦图重心法，单点模型法，相对角法等。

傅健等人通过射线源偏移量和工件上任意一点投影在探测器上的偏移量之间的关系求解各个投影值的位置偏移量，并多次移动射线源，迭代重建CT图像来确定图像重建中心，该方法定位精度不高，而且采用图像信噪比评判图像重建清晰度无法达到预期效果。

张蔚等人依据投影正弦图关于旋转中心对称原理，先对投影正弦图像进行分割，再通过求物体投影在各个方向的投影质心，最后将所有质心位置相加权求均值的方法，此方法优点是降低噪声影响，但是当被扫描物体密度较低时，投影图像分割不准确，导致定位精度较低。

李保磊等人利用Ⅲ代CT扫描、360°旋转方式下投影正弦图相对于旋转中心对称的关系，即正弦图上投影值在旋转中心左侧的面积分与右侧的面积分相等的关系求旋转中心位置，该方法计算精度为像素级，而且受噪声影响大。

LiuT等人运用单点模型法计算旋转中心偏移量，该方法将一根钢丝固定在旋转台上，转台绕旋转中心转360°，运用钢丝投影数据的极值来求解图像重建中心，该模型的缺点是无法避免体积效应引入的误差。

另外一种基于正弦图的相对角法是利用正弦图中旋转角度相差180°的投影值相等特点，以角度相差180°两组投影之差的均值所对应最小的位置作为重建中心，该方法无法避免因探测器具有一定宽度而引入的体积效应，精度不高，并且该方法要求射线源剂量非常稳定

发明内容

本发明的目的就是提供一种基于点扩展函数的工业CT图像重建中心定位方法，它可以确定旋转中心与理想中心的偏移量，从而为修正CT系统旋转中心偏移提供准确的数据。

本发明的目的是通过这样的技术方案实现的，具体步骤如下：

1）将标准待测圆盘工件放置在转台中心，采用Ⅲ代CT扫描方式获取待测圆盘厚度方向中间部位断层投影数据，得到投影正弦图；

2）提取投影正弦图的边缘，计算边缘的质心位置x₀，即图像重建中心初始偏移量；

3）设圆盘旋转中心偏移量x为x₀的±δ领域，即x在[x₀-δ,x₀+δ]范围内，其中0<δ<1，且x的取值为从x₀-δ开始依次递增k像素的序列，再依次运用滤波反投影重建算法重建各个偏移量x对应的CT图像；

4）求出重建后每一幅CT图像，即x每一次取值，中圆盘边缘的点扩展函数PSF的高度h，用一维高斯函数拟合图像重建中心偏移量x与圆盘边缘的点扩展函数PSF的高度h之间的函数关系，计算高斯函数峰值所对应的变量x=b即为所求的图像重建中心偏移量。

进一步，所述待测圆盘质地均匀，采用Ⅲ代CT扫描时扇形射线束平等于圆盘上、下面并垂直于圆盘轴心线，圆盘直径小于等于CT系统扫描视场直径的2/3，且圆盘直径小于等于CT扫描最大可穿透钢厚度的2/3。

进一步，步骤1）中所述的采用Ⅲ代CT扫描方式进行扫描时，旋转台静止，而射线源和探测器组合旋转运动，或旋转台旋转运动，而射线源和探测器组合静止。

进一步，步骤2）中所述的计算边缘的质心位置x₀具体步骤如下：

2-1）采用Canny算法提取投影正弦图的边缘轮廓，其中选取阈值保证所提取的边缘是单像素、连续的轮廓；

2-2）计算第i个投影视下边缘点的重心位置x_i，其中i=1,2,...M，M为投影视角总数；

2-3）计算所有视角下边缘点重心的均值

2-4）计算图像旋转中心初始偏移量其中N为中间探测器位置。

进一步，步骤3）中对不同的图像重建中心偏移量采用FBP算法重建CT图像的总数为1+2δ/k幅，k的取值比系统期望的精度误差低一个数量级。

进一步，步骤4）中根据已重建的1+2δ/k幅CT图像，计算高斯函数峰值所对应的变量x=b即为所求的图像重建中心偏移量的具体步骤如下：

4-1）采用Canny算法提取圆盘CT图像边缘，提取CT图像中圆盘的单像素、封闭的边缘轮廓；

4-2）根据圆盘CT图像边缘，采用最小二乘法求圆盘圆心位置O(x₀,y₀)；

4-3）取两个固定位置A(x₁,y₁)和B(x₂,y₂)，其中点A位于圆盘边缘内部，点B位于圆盘边缘外部；

4-4）求A、B所在圆环半径，即点A、B到O的距离，分别记为r₁、r₂；

4-5）求不同半径r对应的圆环上灰度均值f(r)，r为从r₁到r₂依次以步长1递增的序列；

4-6）以半径r为变量，灰度均值f(r)为因变量作曲线，该曲线即为圆盘CT图像的边缘响应函数ERF；

4-7）对边缘响应函数求一阶导数，即得到PSF曲线，PSF曲线的最大值即为PSF的高度h；

4-8）用一维高斯函数拟合图像重建中心偏移量x与图像PSF高度h之间关系，求出一维高斯函数的三个参数a、b和c；

4-9）高斯函数h(x)峰值对应的变量x=b即为图像重建中心偏移量。

进一步，步骤4-1）包括有以下步骤：

4-1-1）将圆盘CT图像转化为灰度图像；

4-1-2）采用Canny算法提取圆盘图像边缘轮廓，其中选取阈值保证所提取的边缘只包含圆盘边缘，不包含其他噪声边缘，且获取的边缘是封闭的、单像素轮廓；

4-1-3）将边缘位置的行坐标和列坐标依次保存到矩阵X和Y中；

进一步，步骤4-2）所述采用最小二乘法求圆盘圆心位置O(x₀,y₀)的具体步骤如下：

4-2-1）用二元二次函数x²+y²+a₁x+b₁y+c₁=0拟合边缘位置X和Y，其中X和Y分别为常用坐标系中的横坐标和纵坐标；

4-2-2）采用最小二乘拟合求取x²+y²+a₁x+b₁y+c₁=0中参数a₁、b₁和c₁；

4-2-3）计算出圆心位置O(x₀,y₀)为(a₁/2,b₁/2)。

进一步，步骤4-5）的包括以下步骤：

4-5-1）将半径为r₁的圆环上所有像素点对应灰度值的均值作为半径r₁所对应的灰度值；

4-5-2）半径r每增加一个像素，求取一次半径为r圆周上的灰度均值；

4-5-3）当半径r>r₂+1时，停止计算。

由于采用了上述技术方案，本发明具有如下的优点：

本发明所需工件制作简单，操作方便，不受图像环形伪影和硬化伪影的影响，具有高精度、自动定位的优点，并且只需要采样一次数据就能精确定位图像旋转中心，具有很高的实用价值。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书和权利要求书来实现和获得。

附图说明

本发明的附图说明如下。

图1为本发明的流程示意图；

图2为本发明不同半径圆周上的灰度均值示意图；

图3为本发明的图像边缘响应函数曲线；

图4为本发明的图像点扩展函数曲线；

图5为本发明中图像重建中心偏移量与图像点扩展函数高度之间的关系曲线图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

Ⅲ代CT扫描图像重建是基于这样一种假设，扇形射束在旋转投影时，中心射线经过旋转中心，或者说旋转中心在扇形射线束中的位置是已知的，实际CT系统要做到这一点是很困难的，在设备最初调试时要经过大量的试验才能做到；而且设备使用一段时间，因环境和负载变化等因素，这种相对位置也会发生偏移，找准位置费力费时。本方法的目的就是为了实现快速、自动化、更准确地确定旋转中心在扇形射束中相对位置从而确定图像重建中心，实质就是获取实际旋转中心与理想中心的偏移量，从而为修正工件旋转中心偏移提供准确的数据。

一种基于点火器扩展函数的工业CT图像重建中心的方法，如图1所示，包括以下步骤：

1）将标准待测圆盘工件放置在转台中心，采用Ⅲ代CT扫描方式获取待测圆盘厚度方向中间部位断层投影数据，得到投影正弦图；

2）提取投影正弦图的边缘，计算边缘的质心位置x₀，即图像重建中心初始偏移量；

3）设圆盘旋转中心偏移量x为x₀的±δ领域，即x在[x₀-δ,x₀+δ]范围内，其中0<δ<1，且x的取值为从x₀-δ开始依次递增k像素的序列，再依次运用滤波反投影重建算法重建各个偏移量x对应的CT图像；

4）求出重建后每一幅CT图像，即x每一次取值，中圆盘边缘的点扩展函数PSF的高度h，用一维高斯函数拟合图像重建中心偏移量x与圆盘边缘的点扩展函数PSF的高度h之间的函数关系，计算高斯函数峰值所对应的变量x=b即为所求的图像重建中心偏移量。

建立标准试验模型：采用6MeV高能工业CT系统扫描标准圆盘，投影视角数为1024，采样时间为80ms，准直器宽度为0.4mm，CT图像中每个像素代表实际尺寸为0.532mm。

步骤一、将标准圆盘工件放在CT系统转台中心，通过Ⅲ代CT扫描获取圆盘投影正弦图，步骤如下：

1、将标准圆盘工件放置在CT系统转台中心，且尽量使得圆盘中心靠近旋转台中心。

2、采用Ⅲ代CT扫描圆盘厚度方向中间部位，其中CT系统的射线源和探测器保持静止，转台旋转360°，由探测采集与传输子系统得到投影数据即为圆盘投影正弦图。

步骤二、提取投影正弦图的边缘，计算边缘的重心位置x₀,步骤如下：

1、读取投影数据并保存到大小为377×1024的矩阵P中，其中377代表每一投影视角下的投影射线总数，也称为投影探测器个数，1024代表投影视角总数；

2、将投影正弦图P转化为灰度图像，并采用Canny算法提取灰度图像边缘得到P₁；

3、求矩阵P₁第i行边缘点重心位置x_i，i＝1，2，...，1024；

4、通过公式计算P₁的边缘点均值

5、理想重建中心位置为189，图像重建中心的偏移量为x₀=188.9882-189=-0.0118，为了计算方便，x₀取小数点后两位有效数字，即x₀=-0.02；

步骤三、将旋转中心偏移量x设置为x₀的δ邻域[x₀-δ，x₀+δ]内间隔k像素的值，并采用FBP算法逐个重建各个给定偏移量对应的CT图像，步骤如下：

1、取δ=0.5，k=0.01，将图像重建中心偏移量值设置为区间[-0.52，0.48]内从-0.52开始依次递增0.01像素的值；

2、采用基于点驱动的FBP重建算法，依次重建各个偏移量对应的CT图像，为保证精度重建结果以浮点数据格式保存，共重建101幅CT图像。

步骤四、求出重建后每一幅CT图像中圆盘边缘的PSF高度h，用一维高斯函数拟合图像重建中心偏移量与圆盘边缘PSF的高度h之间的函数关系，计算高斯函数峰值所对应的变量x=b即为所求的重建图像中心偏移量，步骤如下：

1、读取图像重建中心偏移量x=-0.02对应的CT图像，采用Canny算法提取图像边缘；

2、根据圆盘图像边缘轮廓的横坐标X和纵坐标Y，采用最小二乘法求取二元二次函数x²+y²+a₁x+b₁y+c₁=0的三个参数a₁、b₁和c₁；计算出圆心位置O(x₀,y₀)为(a₁/2,b₁/2),并求得圆盘的半径为r₀（pixel）；

3、取两个固定位置A(x₁,y₁)和B(x₂,y₂)，其中点A位于圆盘边缘内部，点B位于圆盘边缘外部，且 r₀为圆盘半径，取点示意图如图2所示。

4、按以下方式分别求点A、B所在的圆周半径r₁，r₂，并对r₁，r₂进行取整。

$r_1=\sqrt{{(x_1-x_0)}^2+{(y_1-y_0)}^2}$

$r_2=\sqrt{{(x_2-x_0)}^2+{(y_2-y_0)}^2}$

5、求半径r在r₁～r₂范围内各个圆周上的灰度均值，步骤如下：

1）搜索整幅图像上到圆心距离为r的所有点；

2）统计在半径为r的圆周上的点数n，及其灰度总和sum_I，则半径r对应的灰度均值为f(r)=sum_I/n；

3）计算当前CT图像中不同r所对应的灰度均值f(r)，r∈[r₁,r₂]；

6、以半径r为变量，灰度均值f(r)为因变量，所作曲线即为图像的ERF曲线，如图3所示；

7、边缘响应函数关于半径r求一阶导数，即得到PSF，PSF曲线如图4所示，求PSF曲线最大值即为PSF的高度h；

8、计算101幅CT图像对应的PSF的高度h，表1为不同的旋转中心偏移量x对应的圆盘边缘PSF的高度h；

表1

9、用一维高斯函数拟合表1中偏移量与PSF高度，拟合结果如图5所示，拟合得到三个参数值为a=38.7，b=0.002373，c=1.296，其中拟合误差平方和（SSE）为0.07061，拟合均方根（RMSE）为0.02684，确定系数（R-square）为0.9997；

10、计算得到高斯函数峰值对应的位置x=0.0023727，因此系统图像重建中心偏移量为x≈0.002像素，由于实验中用于数据测试的CT系统其旋转中心已非常逼近理想位置，因此所测试的图像重建中心偏移量应为0，由此可见本发明提供的方法绝对误差小于0.05像素。

为进一步验证本发明中提出方法的有效性，另外采用4MeV高能工业CT系统扫描另一个标准圆盘，投影视角数为1024，采样时间为30ms，准直器宽度为0.2mm，CT图像中每个像素代表实际尺寸为0.186mm。对不同的图像重建中心偏移量分别进行测试，测试结果如表2所示。

表2不同重建中心偏移量测量结果

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (6)

1.一种基于点扩展函数的工业CT图像重建中心定位方法，其特征在于，具体步骤如下：

1)将标准待测圆盘工件放置在转台中心，采用Ⅲ代CT扫描方式获取待测圆盘厚度方向中间部位断层投影数据，得到投影正弦图；

2)提取投影正弦图的边缘，计算边缘的质心位置x₀，即图像重建中心初始偏移量；

3)设图像重建中心偏移量x在x₀的±δ邻域内，即x在[x₀-δ,x₀+δ]范围内，其中0＜δ＜1，且x的取值为从x₀-δ开始依次递增k像素的序列,k的取值比系统期望的精度误差低一个数量级，再依次运用滤波反投影重建算法重建各个偏移量x对应的CT图像；

4)求出重建后每一幅CT图像，即x每一次取值，圆盘边缘的点扩展函数PSF的高度h，用一维高斯函数拟合图像重建中心偏移量x与圆盘边缘的点扩展函数PSF的高度h之间的函数关系，其中a、b和c为一维高斯函数的三个参数，计算高斯函数峰值所对应的变量x＝b即为所求的图像重建中心偏移量；

所述待测圆盘质地均匀，采用Ⅲ代CT扫描时扇形射线束平行于圆盘上、下面并垂直于圆盘轴心线，圆盘直径小于等于CT系统扫描视场直径的2/3，且圆盘直径小于等于CT扫描最大可穿透厚度的2/3；

步骤1)中所述的采用Ⅲ代CT扫描方式进行扫描时，旋转台静止，而射线源和探测器组合旋转运动，或旋转台旋转运动，而射线源和探测器组合静止；

步骤2)中所述的计算边缘的质心位置x₀具体步骤如下：

2-1)采用Canny算法提取投影正弦图的边缘轮廓，其中选取阈值保证所提取的边缘是单像素、连续的轮廓；

2-2)计算第i个投影视角下边缘点的重心位置x_i，其中i＝1,2,...,M，M为投影视角总数；

2-3)计算所有视角下边缘点重心的均值

2-4)计算图像重建中心初始偏移量其中N为中间探测器位置。

2.如权利要求1所述的一种基于点扩展函数的工业CT图像重建中心定位方法，其特征在于，步骤3)中对不同的图像重建中心偏移量采用FBP算法重建CT图像的总数为1+2δ/k幅，k的取值比系统期望的精度误差低一个数量级。

3.如权利要求2所述的一种基于点扩展函数的工业CT图像重建中心定位方法，其特征在于，步骤4)中根据已重建的1+2δ/k幅CT图像，计算高斯函数峰值所对应的变量x＝b即为所求的图像重建中心偏移量的具体步骤如下：

4-1)采用Canny算法提取圆盘CT图像边缘，提取CT图像中圆盘的单像素、封闭的边缘轮廓；

4-2)根据圆盘CT图像边缘，采用最小二乘法求圆盘圆心位置O(x₀,y₀)；

4-3)取两个固定位置A(x₁,y₁)和B(x₂,y₂)，其中点A位于圆盘边缘内部，点B位于圆盘边缘外部；

4-4)求A、B所在圆环半径，即点A、B到O的距离，分别记为r₁、r₂；

4-5)求不同半径r对应的圆环上灰度均值f(r)，r为从r₁到r₂依次以步长1递增的序列；

4-6)以半径r为变量，灰度均值f(r)为因变量作曲线，该曲线即为圆盘CT图像的边缘响应函数ERF；

4-7)对边缘响应函数求一阶导数，即得到PSF曲线，PSF曲线的最大值即为PSF的高度h；

4-8)用一维高斯函数拟合图像重建中心偏移量x与图像PSF高度h之间关系，求出一维高斯函数的三个参数a、b和c；

4-9)高斯函数h(x)峰值对应的变量x＝b即为图像重建中心偏移量。

4.如权利要求3所述的一种基于点扩展函数的工业CT图像重建中心定位方法，其特征在于，步骤4-1)包括有以下步骤：

4-1-1)将圆盘CT图像转化为灰度图像；

4-1-2)采用Canny算法提取圆盘图像边缘轮廓，其中选取阈值保证所提取的边缘只包含圆盘边缘，不包含其他噪声边缘，且获取的边缘是封闭的、单像素轮廓；

4-1-3)将边缘位置的行坐标和列坐标依次保存到矩阵X和Y中。

5.如权利要求4所述的一种基于点扩展函数的工业CT图像重建中心定位方法，其特征在于，步骤4-2)所述采用最小二乘法求圆盘圆心位置O(x₀,y₀)的具体步骤如下：

4-2-1)用二元二次函数x²+y²+a₁x+b₁y+c₁＝0拟合边缘位置X和Y，其中X和Y分别为常用坐标系中的横坐标和纵坐标；

4-2-2)采用最小二乘拟合求取x²+y²+a₁x+b₁y+c₁＝0中参数a₁、b₁和c₁；

4-2-3)计算出圆心位置O(x₀,y₀)为(a₁/2,b₁/2)。

6.如权利要求5所述的一种基于点扩展函数的工业CT图像重建中心定位方法，其特征在于，步骤4-5)包括以下步骤：

4-5-1)将半径为r₁的圆环上所有像素点对应灰度值的均值作为半径r₁所对应的灰度值；

4-5-2)半径r每增加一个像素，求取一次半径为r圆周上的灰度均值；

4-5-3)当半径r＞r₂+1时，停止计算。