CN106408616A - 一种ct成像中投影图背景不一致的校正方法 - Google Patents

Abstract

本发明涉及一种CT成像中投影图背景不一致的校正方法，具体是一种基于图像灰度信息的CT投影图背景互相关配准法，包括如下步骤：1)在采集样品CT投影图像之前，采集亮场图像一张。确定亮场图像的光斑中心；2)计算亮场图像补偿函数；3)亮场图像照度补偿；4)投影图像照度补偿；5)投影图背景与亮场图背景的实际偏移量计算；6)图像平移配准。本发明基于图像灰度信息，利用互相关原理实现对CT投影图的配准，具有不需要预先知道探测器具体偏转参数，操作简单，计算效率高的优点。

Description

一种CT成像中投影图背景不一致的校正方法

技术领域

本发明涉及一种CT成像中投影图背景不一致的校正方法，具体是一种基于图像灰度信息的CT投影图背景互相关配准法。

背景技术

CT成像技术能够快速、无损地探测物体的三维结构，目前已经成为一种重要的材料三维结构表征手段。CT重构过程中，为了避免探测器上的灰尘、缺陷或电子噪声对CT投影图像的影响，经常在CT实验过程中采集亮场图像和暗场图像，通过对投影数据图像扣除背景以提高CT图像精度。高精度的CT图像重构要求在实验过程中保证X射线源、样品台轴心、探测器三者保持相对静止，从而保证探测器在样品转过不同角度时所采集到的投影图具有相同的背景图像。然而在实际成像过程中，由于射线源漂移或机械振动等原因导致的射线源、样品台轴心及探测器三者之间的相对位移很难完全避免，这使得不同角度采集的CT投影图像背景不一致，从而影响CT投影图像的背景扣除以及后续的CT图像重建。

目前，针对CT实验扫描过程中样品台的移动造成CT重建产生的伪影，常用旋转中心测量校正方法校正。常见的旋转中心测量方法有投影边界法、几何关系法、针模法、中心射线法、对称关系法和迭代收敛法等。针对CT实验扫描过程中探测器偏转造成的CT重建图像伪影，常见校正方法有利用原始投影数据正弦图计算校正算法所必需的几何参数、利用含探测器偏转参数的重建公式直接由原始投影数据重建CT图像。但现有方法需要预先知道探测器偏转参数且存在操作复杂和计算量大的技术问题。

发明内容

本发明针对CT扫描重建过程中，由于投影图背景不一致导致的CT重建图像伪影的校正方法存在操作复杂和计算量大的技术问题，提供了一种CT成像中投影图背景不一致的校正方法，这一方法基于互相关原理，不需要预先知道探测器偏转的具体参数，具有操作简便，计算效率高的优点。

为解决上述技术问题，本发明采用的技术方案为：

一种CT成像中投影图背景不一致的校正方法，包括以下步骤：

1)在采集样品CT投影图像之前，采集一张亮场图像，确定亮场图像的光斑中心O(i₀，j₀)；

2)计算亮场图像补偿函数f(r)；

3)亮场图像照度补偿；

4)投影图像照度补偿；

5)投影图背景与亮场图背景的实际偏移量计算；

6)图像平移配准。

所述步骤1)中确定亮场图像的光斑中心O(i₀，j₀)的具体操作步骤为：将亮场图像每点的灰度数值读出，表示为一个灰度数值矩阵{p_i,j}(i＝0,1,2...,i_max；j＝0,1,2,...,j_max)，其中i、j分别表示矩阵{p_i,j}的列坐标和行坐标,与亮场图像上各像素点坐标一一对应；i_max是图像列方向像素点的最大坐标值，j_max是图像行方向像素点的最大坐标；p_i,j表示坐标为(i,j)处的灰度数值；

从亮场图像中光斑边界清晰的图像区域选择第j₁,j₂,...,j_q行，共q行的灰度数据，对每一行灰度数据求其对列坐标i的差分f'_j(i)：

其中Δi取自然数，其数值要在保证光斑边界处的f'_j(i)能明显区别于图像其他位置的前提下取最小值；当差分函数取到最小值f'_j(i)_min和最大值f'_j(i)_max时，对应自变量分别记为和根据公式(2)计算得到光斑在亮场图像第j行的中心坐标

其中表示取数值的整数部分；计算出不同j值时的算术平均值作为亮场图像在行方向的光斑中心坐标i₀，即：

其中,表示取数值的整数部分；

采用同样的办法得到亮场图像列方向光斑中心坐标j₀，则亮场图像光斑中心O点坐标确定为O(i₀，j₀)。

所述步骤2)中亮场图像补偿函数的计算方法为：以光斑中心O(i₀，j₀)为圆心，半径为(r+0.5)作圆，其中r＝[1,2,...,min(i₀,j₀)]，单位为一个像素，min(i₀,j₀)表示取i₀,j₀这两个数值中较小的一个数值。r取不同值时，分别求出圆周上所有像素点灰度值的算术平均值，表示半径为(r+0.5)的圆周上所有像素点灰度值的算术平均值；根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r)，即为亮场图像补偿函数。

所述步骤3)中亮场图像照度补偿的具体操作步骤为：根据亮场图像补偿函数f(r)对亮场图像中所有像素点进行照度补偿：亮场图像灰度数值矩阵{p_i,j}中，坐标为(i,j)的像素点与光源中心O(i₀,j₀)的距离表示为R，则：

将r＝R代入亮场图像补偿函数f(r)中，计算得到亮场图像补偿后在坐标(i,j)处像素点的灰度值a_i,j：

对亮场图像中所有像素点执行公式(5)的计算，得到补偿后的亮场图像灰度数值矩阵{a_i,j}，根据{a_i,j}利用计算机编程可输出补偿后的亮场图像。

所述步骤4)中投影图像照度补偿的具体操作步骤为：

将每张投影图像每点的灰度数值读出，表示为一个灰度数值矩阵(m＝0,1,2...,m_max；n＝0,1,2,...,n_max，k＝1,2,...,k_max)，其中m、n分别表示矩阵的列坐标和行坐标,与投影图像上各像素点坐标一一对应,m_max是图像列方向像素点的最大坐标值，n_max是图像行方向像素点的最大坐标值，k表示投影图的编号，k_max是投影图总的张数；表示第k张投影图上坐标为(m,n)处的灰度数值；确定每张投影图的光斑中心坐标分别表示第k张投影图光斑中心的行坐标和列坐标；利用亮场图像补偿函数f(r)对每张投影图的灰度数值进行补偿，表示补偿后的第k张投影图的灰度数值矩阵，表示照度补偿后的第k张投影图在列坐标为m，行坐标为n处的灰度数值；根据利用计算机编程可输出照度补偿后的投影图像。

所述确定每张投影图的光斑中心坐标的具体操作步骤为：

从第k张投影图像中光斑边界清晰的图像区域选择第n₁,n₂,...,n_z行，共z行的灰度数据，对每一行灰度数据求其对列坐标m的差分

其中Δm取自然数，其数值要在保证光斑边界处的能明显区别于图像其他位置的前提下取最小值；当差分函数取到最小值和最大值时，对应自变量分别记为和根据公式(7)计算得到光斑在第k张投影图像第n行的中心坐标

其中表示取数值的整数部分；计算出不同n值时的算术平均值作为第k张投影图像在行方向的光斑中心坐标即：

其中,表示取数值的整数部分；

采用同样的办法得到第k张投影图像列方向光斑中心坐标则第k张投影图像光斑中心O_k点坐标确定为

所述对每张投影图的灰度数值进行补偿的具体操作步骤为：

根据亮场图像补偿函数f(r)对第k张投影图像中所有像素点进行照度补偿：第k张投影图像灰度数值矩阵中，坐标为(m,n)的像素点与光斑中心的距离表示为D，

则：

将r＝D代入亮场图像补偿函数f(r)中，计算得到第k张投影图像照度补偿后在坐标(m,n)处像素点的灰度值

对第k张图像中所有像素点执行公式(10)的计算，得到补偿后的第k张投影图像灰度数值矩阵根据利用计算机编程可输出补偿后的第k张投影图像。

所述步骤5)中投影图背景与亮场图背景的实际偏移量计算方法为：以照度补偿后的亮场图像作为模板图像，照度补偿后的第k张投影图像作为待配准图像，在模板图像和待配准图像中的相同位置处选取大小相同的区域1和区域2，这两个区域所包含像素点的灰度值分别作为两组离散序列ξ(e,f)和χ^k(e,f)，其中e,f分别表示所选图像区域的列坐标和行坐标，e＝0,1,2,...,e_max，f＝0,1,2,...,f_max，e_max和f_max分别表示所选图像区域中列方向和行方向的最大坐标值，其数值根据所选图像区域大小确定；根据公式(11)计算互相关函数c^k(u,v)的数值，

公式(11)中，u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量，u＝0,1,2，...，u_max,v＝0,1,2,...，v_max,其中u_max与v_max分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量，其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定，代表序列ξ(e,f)的算术平均值，代表序列χ^k(e,f)的算术平均值；当c^k(u,v)取到最大值时，对应的背景偏移量(分别记为)即为亮场图像和第k张投影图像的实际背景偏移量；采用相同的方法将所有投影图像背景配准，从而得到每张投影图像背景与亮场图像背景的实际偏移量。

所述步骤6)中图像平移配准的具体操作步骤为：根据照度补偿后每张投影图背景与亮场图像背景的实际偏移量可以将所有投影图像进行平移变换，实现CT投影图像背景与亮场图像背景的配准；根据背景配准后的第k张投影图灰度数值矩阵利用计算机编程可输出背景配准后的投影图。

所述投影图像平移变换方法如公式(12)所示，

其中，k表示照度补偿后的第k张投影图，表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。

所述步骤1)中,所选择的总行数q可以根据图像情况选择不同数目，选择的行数越多，计算得到亮场图像在行方向的光斑中心位置越精确。在确定列方向光斑中心时同样如此。

步骤1)中,在采集样品CT投影图像之前，采集亮场图像一张。从亮场图像中光斑边界清晰的图像区域选择第j₁,j₂,...,j_q行，共q行的灰度数据，对每一行灰度数据求其对列坐标i的差分f'_j(i)。亮场图像采集方法以及确定亮场图像中光斑边界清晰的图像区域为本领域常规方法。

步骤2)中,在根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r)。多项式拟合法为本领域常规使用方法。

步骤4)中,从第k张投影图像中光斑边界清晰的图像区域选择第n₁,n₂,...,n_z行，共z行的灰度数据，对每一行灰度数据求其对列坐标m的差分所选择的总行数z可以根据图像情况选择不同数目。选择的行数越多，计算得到投影图像在行方向的光斑中心位置越精确。在确定列方向光斑中心时同样如此。确定投影图像中光斑边界清晰的图像区域是本领域常规方法。

步骤5)中,在其中u_max与v_max分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量，其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定。确定u_max与v_max的数值是本领域常规使用方法。

本发明中提到的将图像灰度数值读出表示为一个灰度数值矩阵、将灰度数值矩阵输出为图像均是利用计算机编程实现，是本领域常规使用手段。本发明中提到的对灰度数值矩阵中全部或部分坐标处的灰度值根据不同公式进行计算都是利用计算机编程实现，属于本领域常规使用手段。

本发明采用以上技术方案，与背景技术相比，具有以下优点：

1、不需要预先知道探测器偏转参数；

2、操作简单，计算量小。

附图说明

图1为亮场图像；

图2为图1中横线位置处灰度值；

图3为照度补偿后亮场图像；

图4为图3中横线位置处灰度值分布；

图5为第一张投影图；

图6为图5中横线位置处的灰度值分布；

图7为照度补偿后的第一张投影图；

图8为图7中横线位置处灰度值分布；

图9为背景校准后的第一张投影图；

图10为利用未进行背景配准的CT投影图重建得到的某张CT切片；

图11为利用背景配准后的CT投影图重建得到的某张CT切片。

具体实施方式

本发明所述的CT投影图像背景修正方法可修正不同CT设备获取的不同样品的CT投影图像背景。本实施例所用投影数据是在三英精密仪器有限公司制造的实验室CT机上采集到的。为了验证互相关配准法的有效性，数据采集过程中人为引入探测器抖动。样品为直径2.5mm的圆柱状人造砂岩与直径2.85mm的圆柱状纯铝样品粘合在一起。亮场图像、暗场图像、样品CT投影图的像素均为2048×2048，本次实验共采集投影图像901张。该实施例只是为了进一步阐述本发明，但并不限制本发明所保护的范围。

一种CT成像中投影图背景不一致的校正方法，包括以下步骤：

1)在采集样品CT投影图像之前，采集一张亮场图像，确定亮场图像的光斑中心O(i₀，j₀)：将亮场图像每点的灰度数值读出，表示为一个灰度数值矩阵{p_i,j}(i＝0,1,2...,i_max；j＝0,1,2,...,j_max)，其中i、j分别表示矩阵{p_i,j}的列坐标和行坐标,与亮场图像上各像素点坐标一一对应；i_max是图像列方向像素点的最大坐标值，j_max是图像行方向像素点的最大坐标；p_i,j表示坐标为(i,j)处的灰度数值；

从亮场图像中光斑边界清晰的图像区域选择第j₁,j₂,...,j_q行，共q行的灰度数据，对每一行灰度数据求其对列坐标i的差分f_j'(i)：

其中Δi取自然数，其数值要在保证光斑边界处的f_j'(i)能明显区别于图像其他位置的前提下取最小值；当差分函数取到最小值f_j'(i)_min和最大值f_j'(i)_max时，对应自变量分别记为和根据公式(2)计算得到光斑在亮场图像第j行的中心坐标

其中表示取数值的整数部分；计算出不同j值时的算术平均值作为亮场图像在行方向的光斑中心坐标i₀，即：

其中,表示取数值的整数部分；

采用同样的办法得到亮场图像列方向光斑中心坐标j₀，则亮场图像光斑中心O点坐标确定为O(i₀，j₀)。

这一步中,所选择的总行数q可以根据图像情况选择不同数目，选择的行数越多，计算得到亮场图像在行方向的光斑中心位置越精确。在确定列方向光斑中心时同样如此。

图1所示为CT实验中采集的亮场图像，图像像素为2048×2048，将图像中所有点的灰度值读出可以表示为一个灰度数值矩阵{p_i,j}(i＝0,1,2...,2047；j＝0,1,2,...,2047)。图2为图1中横线位置处灰度值，由图2可看到亮场图像灰度值沿图1中横线方向至左向右明显减小，灰度值分布不均匀。本文选取图像中光斑边界清晰的第100-150行(图中a点与b点之间)共51行的灰度数值数据，对每一行灰度数据求其对列坐标i的差分f'_j(i)：

以j＝100(即第100行)为例，若Δi＝1，差分涨落参差不齐，很难找到f'₁₀₀(i)的极值。若Δi＝2，除光斑边界处发生突变外，其他区域f'₁₀₀(i)的涨落均匀，即选Δi＝2较为理想。根据公式(1)计算得到，当差分函数取到最小值f'₁₀₀(i)_min＝-200.5时，对应自变量为光斑在第100行的左边界；当差分函数取到最大值f'₁₀₀(i)_max＝236.75时，对应自变量为光斑在第100行的右边界。根据公式(2)计算得到光斑在亮场图像第100行的中心坐标利用同样的方法计算出图像中第100行到150行中所有行的中心坐标，根据公式(3)计算得到亮场图像在行方向的光斑中心坐标为1018。采用同样的方法，计算得到亮场图像在列方向的光斑中心坐标为1017，这样确定亮场图像光斑中心O点坐标(1018，1007)。

2)计算亮场图像补偿函数f(r)：以光斑中心O(i₀，j₀)为圆心，半径为(r+0.5)作圆，其中r＝[1,2,...,min(i₀,j₀)]，单位为一个像素，min(i₀,j₀)表示取i₀,j₀这两个数值中较小的一个数值。r取不同值时，分别求出圆周上所有像素点灰度值的算术平均值，表示半径为(r+0.5)的圆周上所有像素点灰度值的算术平均值；根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r)，即为亮场图像补偿函数。

以图1所示亮场图像中坐标为(1018，1007)处的像素点为圆心，半径分别取r＝1.5,2.5，...,1007.5，在图像上做圆，并分别计算出每个半径对应圆周上像素点的平均灰度值利用多项式拟合法得到关于r的函数：f(r)＝4.90051r³-0.03556r²+0.0000118r

3)亮场图像照度补偿。根据亮场图像补偿函数f(r)对亮场图像中所有像素点进行照度补偿：亮场图像灰度数值矩阵{p_i,j}中，坐标为(i,j)的像素点与光源中心O(i₀,j₀)的距离表示为R，则：

将r＝R代入亮场图像补偿函数f(r)中，计算得到亮场图像补偿后在坐标(i,j)处像素点的灰度值a_i,j：

对亮场图像中所有像素点执行公式(5)的计算，得到补偿后的亮场图像灰度数值矩阵{a_i,j}，根据{a_i,j}利用计算机编程可输出补偿后的亮场图像。

根据公式(4)计算出图1中除圆心外，其他所有点到圆心O的距离R，将R的数值代入函数f(r)＝4.90051r³-0.03556r²+0.0000118r中，计算出f(R)后，代入公式(5)可计算出照度补偿后(i,j)点处的灰度值a_i,j。对图1中所有点进行上述操作后输出补偿后的亮场图像，如图3所示。图4为图3中横线位置处灰度值分布。通过对比图1图2与图3图4可看出照度补偿后亮场图像灰度分布明显变得较为均匀。

4)投影图像照度补偿：将每张投影图像每点的灰度数值读出，表示为一个灰度数值矩阵(m＝0,1,2...,m_max；n＝0,1,2,...,n_max，k＝1,2,...,k_max)，其中m、n分别表示矩阵的列坐标和行坐标,与投影图像上各像素点坐标一一对应,m_max是图像列方向像素点的最大坐标值，n_max是图像行方向像素点的最大坐标值，k表示投影图的编号，k_max是投影图总的张数；表示第k张投影图上坐标为(m,n)处的灰度数值；确定每张投影图的光斑中心坐标分别表示第k张投影图光斑中心的行坐标和列坐标；利用亮场图像补偿函数f(r)对每张投影图的灰度数值进行补偿，表示补偿后的第k张投影图的灰度数值矩阵，表示照度补偿后的第k张投影图在列坐标为m，行坐标为n处的灰度数值；根据利用计算机编程可输出照度补偿后的投影图像。

确定每张投影图的光斑中心坐标的具体操作步骤为：从第k张投影图像中光斑边界清晰的图像区域选择第n₁,n₂,...,n_z行，共z行的灰度数据，对每一行灰度数据求其对列坐标m的差分

其中Δm取自然数，其数值要在保证光斑边界处的能明显区别于图像其他位置的前提下取最小值；当差分函数取到最小值和最大值时，对应自变量分别记为和根据公式(7)计算得到光斑在投影图像第n行的中心坐标

其中表示取数值的整数部分；计算出不同n值时的算术平均值作为第k张投影图像在行方向的光斑中心坐标即：

其中,表示取数值的整数部分；

采用同样的办法得到第k张投影图像列方向光斑中心坐标则第k张投影图像光斑中心O_k点坐标确定为

对每张投影图的照度进行补偿的具体操作步骤为：根据亮场图像补偿函数f(r)对第k张投影图像中所有像素点进行照度补偿：第k张投影图像灰度数值矩阵中，坐标为(m,n)的像素点与光斑中心的距离表示为D，则：

将r＝D代入亮场图像补偿函数f(r)中，计算得到第k张投影图像照度补偿后在坐标(m,n)处像素点的灰度值

对第k张投影图像中所有像素点执行公式(10)的计算，得到补偿后的第k张投影图像灰度数值矩阵根据利用计算机编程可输出照度补偿后的第k张投影图像。

以CT实验中的第一张投影图为例，如图5所示，图6为图5中横线位置处灰度数值分布，由图6可看出图5中沿横线方向灰度数值至左向右明显减小，分布不均匀。将图5中所有点灰度数值读出可表示为一个灰度数值矩阵(m＝0,1,2...,2047；n＝0,1,2,...,2047)，从第1张投影图像中光斑边界清晰的图像区域选择第100,101,...,150行(图中c点和d点之间)，共51行的灰度数据，对每一行灰度数据求其对列坐标m的差分

以第100行的数据为例，当Δm＝1时，计算得到第100行的差分涨落参差不齐，很难找到的极值。当Δm＝2时，除光斑边界处发生突变外，其他区域的涨落均匀，即选Δm＝2较为理想。当差分函数取到最小值时，对应自变量为光斑左边界；当差分函数取到最大值时，对应自变量为光斑右边界。根据公式(7)计算得到光斑在投影图像第100行的中心坐标利用同样的方法计算出第一张投影图像中第100行到150行中所有行的中心坐标，根据公式(8)计算得到第一张投影图像在行方向的光斑中心坐标为992。采用同样的方法，计算得到第一张投影图像在列方向的光斑中心坐标为1006，这样确定第一张投影图像光斑中心点坐标(992，1006)。

以第一张投影图为例说明投影图像的照度补偿，根据公式(9)计算出图5中除圆心外，其他所有点到圆心O的距离D，将D的数值代入函数f(r)＝4.90051r³-0.03556r²+0.0000118r中，计算出f(D)后，代入公式(10)可计算出照度补偿后(m,n)点处的灰度值对图5中所有点进行上述操作后输出补偿后的第一张投影图像，如图7所示。图8为图7中横线位置处灰度值分布。通过对比图4图5与图7图8可看出经过照度补偿后的第一张投影图像灰度数值分布变得较为均匀。

5)投影图背景与亮场图背景的实际偏移量计算：以照度补偿后的亮场图像作为模板图像，照度补偿后的第k张投影图像作为待配准图像，在模板图像和待配准图像中的相同位置处选取大小相同的区域1和区域2，这两个区域所包含像素点的灰度值分别作为两组离散序列ξ(e,f)和χ^k(e,f)，其中e,f分别表示所选图像区域的列坐标和行坐标，e＝0,1,2,...,e_max，f＝0,1,2,...,f_max，e_max和f_max分别表示所选图像区域中列方向和行方向的最大坐标值，其数值根据所选图像区域大小确定；根据公式(11)计算互相关函数c^k(u,v)的数值，

公式(11)中，u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量，u＝0,1,2，...，u_max,v＝0,1,2,...，v_max,其中u_max与v_max分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量，其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定，代表序列ξ(e,f)的算术平均值，代表序列χ^k(e,f)的算术平均值；当c^k(u,v)取到最大值时，对应的背景偏移量(分别记为)即为亮场图像和第k张投影图像的实际背景偏移量；采用相同的方法将所有投影图像背景配准，从而得到每张投影图像背景与亮场图像背景的实际偏移量。

以照度补偿后的亮场图(图3)为模板，照度补偿后的投影图(图7)为待配准图像，在两张图中相同位置处选取大小相同的区域(图3中的区域1和图7中的区域2)。把这两个区域中的灰度数值读出形成两组离散序列ξ(e,f)和χ^k(e,f)，其中e,f分别表示所选图像区域的列坐标和行坐标，e＝0,1,2,...,34，f＝0,1,2,...,691，k＝1,根据公式(11)计算互相关函数c¹(u,v)的数值。

公式(11)中，u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量，其中u＝0,1,2，...，u_max,v＝0,1,2,...，v_max,根据本实施例中X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定投影图与亮场图像在列方向和行方向可能的最大背景偏移量u_max＝100，v_max＝100。代表序列ξ(e,f)的算术平均值，代表序列χ¹(e,f)的算术平均值。

通过计算发现，当c¹(u,v)取到最大值0.988时，对应的即为亮场图像和第1张投影图像的实际背景偏移量。

6)图像平移配准：根据照度补偿后每张投影图背景与亮场图像背景的实际偏移量可以将所有投影图像进行平移变换，实现CT投影图像背景与亮场图像背景的配准；根据背景配准后的第k张投影图灰度数值矩阵利用计算机编程可输出背景配准后的投影图。投影图像平移变换方法如公式(12)所示，

其中，k表示照度补偿后的第k张投影图，表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。

以照度补偿后的第一张投影图像为例，根据亮场图像和第1张投影图像的实际背景偏移量按照公式(12)对照度补偿后的第一张投影图进行平移变换，因为在第一张投影图中总是有所以平移变换公式(12)的具体形式变为:

其中，表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。对投影图进行平移变换后，灰度数值矩阵即为背景校准后的第1张投影图灰度数值矩阵，利用计算机编程输出该图像如图9所示。图10和图11分别为利用未进行背景配准的CT投影图和背景配准后的CT投影图重建得到样品相同位置处CT切片。对比图10和图11可看出，对投影图像背景配准后明显提高了CT切片重构质量。

Claims (10)

1.一种CT成像中投影图背景不一致的校正方法，其特征在于：包括以下步骤：

1)在采集样品CT投影图像之前，采集一张亮场图像，确定亮场图像的光斑中心O(i₀，j₀)；

2)计算亮场图像补偿函数f(r)；

3)亮场图像照度补偿；

4)投影图像照度补偿；

5)投影图背景与亮场图背景的实际偏移量计算；

6)图像平移配准。

2.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述步骤1)中确定亮场图像的光斑中心O(i₀，j₀)的具体操作步骤为：将亮场图像每点的灰度数值读出，表示为一个灰度数值矩阵{p_i,j}(i＝0,1,2...,i_max；j＝0,1,2,...,j_max)，其中i、j分别表示矩阵{p_i,j}的列坐标和行坐标,与亮场图像上各像素点坐标一一对应；i_max是图像列方向像素点的最大坐标值，j_max是图像行方向像素点的最大坐标；p_i,j表示坐标为(i,j)处的灰度数值；

从亮场图像中光斑边界清晰的图像区域选择第j₁,j₂,...,j_q行，共q行的灰度数据，对每一行灰度数据求其对列坐标i的差分f′_j(i)：

$f_j^{\′}\left(i\right)=\frac{p_{j,i+\mathrm{\Δ}i}-p_{j,i-\mathrm{\Δ}i}}{2\mathrm{\Δ}i},(i=\mathrm{\Δ}i,\mathrm{\Δ}i+1,\mathrm{\Δ}i+\mathrm{2...},i_{max}-\mathrm{\Δ}i;j=j_1,j_2,\mathrm{...},j_q)---\left(1\right)$

其中Δi取自然数，其数值要在保证光斑边界处的f′_j(i)能明显区别于图像其他位置的前提下取最小值；当差分函数取到最小值f′_j(i)_min和最大值f′_j(i)_max时，对应自变量分别记为和根据公式(2)计算得到光斑在亮场图像第j行的中心坐标

$i_0^j=floor(\frac{i_1^j+i_2^j}{2}+0.5)---\left(2\right)$

其中表示取数值的整数部分；计算出不同j值时的算术平均值作为亮场图像在行方向的光斑中心坐标i₀，即：

$i_0=floor\left(\frac{1}{q}\underset{j=j_1}{\overset{j_q}{\Σ}}{i}_0^{\left(j\right)}\right)---\left(3\right)$

其中,表示取数值的整数部分；

采用同样的办法得到亮场图像列方向光斑中心坐标j₀，则亮场图像光斑中心O点坐标确定为O(i₀，j₀)。

3.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述步骤2)中亮场图像补偿函数的计算方法为：以光斑中心O(i₀，j₀)为圆心，半径为(r+0.5)作圆，其中r＝[1,2,...,min(i₀,j₀)]，单位为一个像素，min(i₀,j₀)表示取i₀,j₀这两个数值中较小的一个数值。r取不同值时，分别求出圆周上所有像素点灰度值的算术平均值，表示半径为(r+0.5)的圆周上所有像素点灰度值的算术平均值；根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r)，即为亮场图像补偿函数。

4.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述步骤3)中亮场图像照度补偿的具体操作步骤为：根据亮场图像补偿函数f(r)对亮场图像中所有像素点进行照度补偿：亮场图像灰度数值矩阵{p_i,j}中，坐标为(i,j)的像素点与光源中心O(i₀,j₀)的距离表示为R，则：

$R={({\left(i-i_0\right)}^2+{\left(j-j_0\right)}^2)}^{\frac{1}{2}}---\left(4\right)$

将r＝R代入亮场图像补偿函数f(r)中，计算得到亮场图像补偿后在坐标(i,j)处像素点的灰度值a_i,j：

$a_{i,j}=\frac{p_{i,j}}{f\left(R\right)}\×40000,(i=0,1,\mathrm{2...},i_{max};j=0,1,2,\mathrm{...},j_{max})---\left(5\right)$

对亮场图像中所有像素点执行公式(5)的计算，得到补偿后的亮场图像灰度数值矩阵{a_i,j}，根据{a_i,j}利用计算机编程可输出补偿后的亮场图像。

5.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述步骤4)中投影图像照度补偿的具体操作步骤为：

将每张投影图像每点的灰度数值读出，表示为一个灰度数值矩阵(m＝0,1,2...,m_max；n＝0,1,2,...,n_max，k＝1,2,...,k_max)，其中m、n分别表示矩阵的列坐标和行坐标,与投影图像上各像素点坐标一一对应,m_max是图像列方向像素点的最大坐标值，n_max是图像行方向像素点的最大坐标值，k表示投影图的编号，k_max是投影图总的张数；表示第k张投影图上坐标为(m,n)处的灰度数值；确定每张投影图的光斑中心坐标 分别表示第k张投影图光斑中心的行坐标和列坐标；利用亮场图像补偿函数f(r)对每张投影图的灰度数值进行补偿，表示补偿后的第k张投影图的灰度数值矩阵，表示照度补偿后的第k张投影图在列坐标为m，行坐标为n处的灰度数值；根据利用计算机编程可输出照度补偿后的投影图像。

6.根据权利要求5所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述确定每张投影图的光斑中心坐标的具体操作步骤为：

从第k张投影图像中光斑边界清晰的图像区域选择第n₁,n₂,...,n_z行，共z行的灰度数据，对每一行灰度数据求其对列坐标m的差分

${f^{\′}}_n^k\left(m\right)=\frac{T_{n,m+\mathrm{\Δ}m}^k-T_{n,m-\mathrm{\Δ}m}^k}{2\mathrm{\Δ}m},(m=\mathrm{\Δ}m,\mathrm{\Δ}m+1,\mathrm{\Δ}m+\mathrm{2...},m_{\max }-\mathrm{\Δ}m;n=n_1,n_2,\mathrm{...},n_z)---\left(6\right)$

其中Δm取自然数，其数值要在保证光斑边界处的能明显区别于图像其他位置的前提下取最小值；当差分函数取到最小值和最大值时，对应自变量分别记为和根据公式(7)计算得到光斑在第k张投影图像第n行的中心坐标

$m_0^{k,n}=floor(\frac{m_1^{k,n}+m_2^{k,n}}{2}+0.5)---\left(7\right)$

其中表示取数值的整数部分；计算出不同n值时的算术平均值作为第k张投影图像在行方向的光斑中心坐标即：

$m_0^k=floor\left(\frac{1}{z}\underset{n=n_1}{\overset{n_z}{\Σ}}{m}_0^{k,n}\right)---\left(8\right)$

其中,表示取数值的整数部分；

采用同样的办法得到第k张投影图像列方向光斑中心坐标则第k张投影图像光斑中心O_k点坐标确定为

7.根据权利要求5所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述对每张投影图的灰度数值进行补偿的具体操作步骤为：

根据亮场图像补偿函数f(r)对第k张投影图像中所有像素点进行照度补偿：第k张投影图像灰度数值矩阵中，坐标为(m,n)的像素点与光斑中心的距离表示为D，则：

$D={({\left(m-m_0^k\right)}^2+{\left(n-n_0^k\right)}^2)}^{\frac{1}{2}}---\left(9\right)$

将r＝D代入亮场图像补偿函数f(r)中，计算得到第k张投影图像照度补偿后在坐标(m,n)处像素点的灰度值

$b_{m,n}^k=\frac{T_{m,n}}{f\left(D\right)}\×40000,(m=0,1,2,\mathrm{...},m_{max};n=0,1,2,\mathrm{...},n_{max})---\left(10\right)$

对第k张图像中所有像素点执行公式(10)的计算，得到补偿后的第k张投影图像灰度数值矩阵根据利用计算机编程可输出补偿后的第k张投影图像。

8.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述步骤5)中投影图背景与亮场图背景的实际偏移量计算方法为：以照度补偿后的亮场图像作为模板图像，照度补偿后的第k张投影图像作为待配准图像，在模板图像和待配准图像中的相同位置处选取大小相同的区域1和区域2，这两个区域所包含像素点的灰度值分别作为两组离散序列ξ(e,f)和χ^k(e,f)，其中e,f分别表示所选图像区域的列坐标和行坐标，e＝0,1,2,...,e_max，f＝0,1,2,...,f_max，e_max和f_max分别表示所选图像区域中列方向和行方向的最大坐标值，其数值根据所选图像区域大小确定；根据公式(11)计算互相关函数c^k(u,v)的数值，

$c^k(u,v)=\frac{\underset{e=0,f=0}{\overset{e=e_{\max },f=f_{\max }}{\Σ}}\left(\right(\mathrm{\ξ}\left(e,f\right)-\stackrel{\‾}{\mathrm{\ξ}}\left)\right({\mathrm{\χ}}^k\left(e-u,f-v\right)-\stackrel{\‾}{{\mathrm{\χ}}^k}\left)\right)}{{\left(\underset{e=0,f=0}{\overset{e=e_{\max },f=f_{\max }}{\Σ}}{\left(\mathrm{\ξ}\left(e,f\right)-\stackrel{\‾}{\mathrm{\ξ}}\right)}^2\right)}^{\frac{1}{2}}{\left(\underset{e=0,f=0}{\overset{e=e_{\max },f=f_{\max }}{\Σ}}{\left({\mathrm{\χ}}^k\left(e-u,f-v\right)-\stackrel{\‾}{{\mathrm{\χ}}^k}\right)}^2\right)}^{\frac{1}{2}}}---\left(11\right)$

公式(11)中，u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量，u＝0,1,2，...，u_max,v＝0,1,2,...，v_max,其中u_max与v_max分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量，其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定，代表序列ξ(e,f)的算术平均值，代表序列χ^k(e,f)的算术平均值；当c^k(u,v)取到最大值时，对应的背景偏移量(分别记为)即为亮场图像和第k张投影图像的实际背景偏移量；采用相同的方法将所有投影图像背景配准，从而得到每张投影图像背景与亮场图像背景的实际偏移量。

9.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述步骤6)中图像平移配准的具体操作步骤为：根据照度补偿后每张投影图背景与亮场图像背景的实际偏移量可以将所有投影图像进行平移变换，实现CT投影图像背景与亮场图像背景的配准；根据背景配准后的第k张投影图灰度数值矩阵利用计算机编程可输出背景配准后的投影图。

10.根据权利要求9所述的一种CT成像中投影图背景不一致的校正方法，其特征在于：所述投影图像平移变换方法如公式(12)所示，

$\left\{\begin{array}{cc}{b^{\&prime;}}_{m,n}^k=b_{m-u_0^k,n-v_0^k}^k& m\&GreaterEqual;u_0^k,n\&GreaterEqual;v_0^k\\ {b^{\&prime;}}_{m,n}^k=b_{m,n-v_0^k}^k& m<u_0^k,n\&GreaterEqual;v_0^k\\ {b^{\&prime;}}_{m,n}^k=b_{m-u_0^k,n}^k& m\&GreaterEqual;u_0^k,n<v_0^k\end{array}\right.---\left(12\right)$

其中，k表示照度补偿后的第k张投影图，表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。