CN101718719A - 一种连续扫描三维锥束工业ct角度增量确定方法 - Google Patents

Abstract

本发明属于X射线计算机层析成像(CT)技术领域，具体为一种连续扫描三维锥束工业CT角度增量确定方法，步骤为：(1)进行单圆轨道锥束CT扫描，获得一组二维投影图像；(2)对步骤(1)中二维投影图像进行对数解调，获得一组二维线积分图像；(3)计算步骤(2)中第一幅二维线积分图像与步骤(2)中其他二维线积分图像的相关系数，形成一个一维相关系数数组；(4)搜索步骤(3)获得的一维相关系数数组中的最大值，记录下该最大值在该数组中的序号；(5)将步骤(4)获得的最大值对应的序号减1，并除2π，结果即为连续扫描三维锥束工业CT的角度增量。本发明实现过程简单、高效，精度高，不需要特殊的硬件和额外的锥束扫描投影数据。

Description

一种连续扫描三维锥束工业CT角度增量确定方法

技术领域

本发明涉及一种连续扫描三维锥束工业CT角度增量确定方法，属于X射线计算机层析成像(CT)技术领域。

背景技术

在X射线CT系统中，X射线源发出X射线，从不同角度穿过被检测物体的某一区域，放置于射线源对面的探测器在相应角度接受，然后根据各角度射线不同程度的衰减，利用一定的重建算法和计算机进行运算，重建出物体被扫描区域的射线线衰减系数分布映射图像，从而实现由投影重建图像，无损地再现物体在该区域内的介质密度、成分和结构形态等特征。

因为较高的成像效率，基于面阵探测器的X射线三维锥束工业CT成像技术在航空、航天、核工业等领域得到了越来越广泛的应用。通常，这类系统采用单圆轨道锥束扫描方式实施CT扫描。在这种扫描方式下，X射线源和探测器是静止的，而物体放置于转台，在360°范围内步进旋转；探测器在每个扫描角度下采集透过物体的X射线投影信号，形成一组二维投影图像；最后，基于这组二维投影图像和FDK类型的重建算法重建物体三维CT图像。这类系统的成像速度主要由扫描速度和重建速度决定。目前，由于图形处理单元(GPU)及其相应并行运算架构CUDA的应用，锥束CT的重建速度已显著提高，大大高于扫描速度。因此，为进一步提高成像效率，必须减小扫描时间。减小扫描时间的一种方法是用连续扫描方式取代步进扫描方式。

在实际应用中，与步进扫描方式相比，连续扫描方式存在一个突出问题：由于难以做到探测器和转台的工作时序完全同步，CT扫描角度增量不能准确确定。U.Kumar等，A statistical correction method for minimization of systemic artefact in acontinuous-rotate X-ray based industrial CT system，Nuclear Instruments andMethods in Physics Research A，Vol.515，pp.829-839，2003，描述了一种用于连续扫描二维工业CT角度增量确定方法。这种方法首先计算二维工业CT投影正弦图第一行数据与其他各行数据差值的标准差，然后以最小标准差对应的扫描位置为完成360°扫描的位置，进而确定扫描角度增量。但是，这种方法仅适用于基于线阵探测器的二维工业CT系统，难以为基于面阵探测器的连续扫描三维锥束工业CT系统采用。

目前，尚未发现有用于基于面阵探测器的连续扫描三维锥束工业CT系统的角度增量确定方法。

发明内容

本发明要解决的技术问题是：针对基于面阵探测器的连续扫描三维锥束工业CT系统角度增量确定问题，提供一种基于投影图像相关性的角度增量确定方法。该方法实现过程简单、高效，精度高，不需要特殊的硬件和额外的锥束扫描投影数据。

本发明采用的技术方案：连续扫描三维锥束工业CT角度增量确定方法，其特征在于包括下列步骤：

(1)进行单圆轨道锥束CT扫描，获得一组二维投影图像；

(2)对步骤(1)中所述二维投影图像进行对数解调，获得一组二维线积分图像；

(3)计算步骤(2)中第一幅二维线积分图像与步骤(2)中其他二维线积分图像的相关系数，形成一个一维相关系数数组；

(4)搜索步骤(3)获得的一维相关系数数组中的最大值，记录下该最大值在该数组中的序号；

(5)将步骤(4)获得的最大值对应的序号减1，并除2π，结果即为连续扫描三维工业CT的角度增量。

本发明基于面阵探测器的连续扫描三维锥束工业CT扫描原理如图1：扫描物体放置于检台上匀速连续旋转；旋转过程中，由面阵探测器以固定采样速度连续采集透射过物体的射线投影，获得二维投影图像；当检台旋转380度时，探测器停止采样，检台和射线源同时停止，扫描结束；根据扫描形成的投影数据和角度增量，利用对数解调算法和CT重建算法进行重建，即可获得物体三维层析图像。

在如图1所示的连续扫描方式下，探测器获取的二维投影图像可用I(i，m，n)表示，其中，i表示扫描角度，(m，n)表示某一探测通道在面阵探测器上的位置。根据这些二维投影图像，可以用对数解调公式(1)获得相应的二维线积分图像P(i，m，n)。

$P(i,m,n)=\ln \left(\frac{\mathrm{mean}\left(I(i,1:\mathrm{10,1}:K)\right)}{I(i,m,n)}\right)---\left(1\right)$

其中，P(i，m，n)表示第i个扫描角度对应的二维线积分图像；ln表示自然对数运算；mean表示二维均值运算；1∶10表示从1取到10，1∶K表示从1取到K，K为二维线积分图像的高度，(1∶10，1∶K)定义了二维图像上的一个宽度为10高度为K的区域。例如，设第i个扫描角度对应的二维投影图像I(i，m，n)的宽度为900，高度为1000，则其二维线积分图像求取公式为 $P(i,m,n)=\ln \left(\frac{\mathrm{mean}\left(I(i,1:\mathrm{10,1}:1000)\right)}{I(i,m,n)}\right).$

利用FDK类型的CT重建算法进行重建，除了要利用上述二维线积分图像，还必须计算出正确的连续角度增量。即使是轻微的角度增量错误，在CT图像上也会导致明显的伪影和模糊。为说明这个问题，我们利用经典三维Shepp-Logan模型进行了计算机仿真研究。图2(a)显示了一组没有角度误差的CT重建图像，图2(b)显示了角度误差为0.0028°时的CT重建图像，图2(c)显示了角度误差为0.0056°时的CT重建图像。图3(a)是图2(a)与图2(b)的差值图像，图3(b)是图2(a)与图2(c)的差值图像。图3(c)显示了图2(a)、(b)和(c)中线条对应部分的灰度曲线。从这些结果可以看出，轻微的角度增量错误也会在CT图像上导致明显的伪影和模糊。

为此，本发明提出一种基于面阵探测器的连续扫描三维锥束工业CT系统角度增量确定方法。该方法利用了CT扫描的如下性质：根据CT扫描原理，同一物体在0°和360°扫描角度下，其二维线积分图像相同。它首先计算0°扫描角的二维线积分图像与其他扫描角度下的二维线积分图像的相关系数，然后将具有最大相关系数的扫描角度视为360°扫描角位置。当扫描角度360°已知后，就可根据0°和360°下的二维线积分图像的序号利用公式(2)计算出360°内探测器采集的图像的个数。

C(i)＝corr(P(1，m，n)，P(i，m，n))，1≤i≤N_p

N＝i|_{C(i)＝max(C)}    (2)

公式(2)中，C(i)表示第1个二维线积分图像与第i个二维线积分图像的相关系数；corr表示二维相关运算；P(1，m，n)表示第1个扫描角度对应的二维线积分图像；P(i，m，n)表示第i个扫描角度对应的二维线积分图像；N_p表示物体在380度内采集的所有图像的个数；N表示与第1个二维线积分图像具有最大相关性的二维线积分图像对应的扫描角度序号，也表示360°内探测器采集的图像的个数；C表示第1个二维线积分图像与其他二维线积分图像的相关系数；max表示一维数组取最大值运算。

当360°内探测器采集的图像的个数N已知后，即可利用公式(3)计算出连续扫描角度增量。这里，A表示角度增量。

A＝2π/(N-1)    (3)

本发明与现有技术相比的优点如下：

(1)本发明解决了连续扫描三维锥束CT角度增量确定问题；

(2)本发明对整个二维图像进行相关运算，可降低数据中统计噪声对结果的影响，故精度高；

(3)本发明实现过程简单、高效，不需要特殊的硬件和额外的锥束扫描投影数据，易于工程实现。

附图说明

图1为基于面阵探测器的连续扫描三维锥束工业CT扫描结构图；

图2(a)为一组没有角度误差的CT重建图像；

图2(b)为角度误差为0.0028°时的CT重建图像；

图2(c)为角度误差为0.0056°时的CT重建图像；

图3(a)是图2(a)与图2(b)的差值图像；

图3(b)是图2(a)与图2(c)的差值图像；

图3(c)是图2(a)、(b)和(c)中线条对应部分的灰度曲线；

图4是相关系数曲线图；

图5(a)是第1幅二维投影图像；

图5(b)是第902幅二维投影图像；

图5(c)是第960幅二维投影图像；

图5(d)是第1幅与第902幅二维投影图像的差值图像；

图5(e)是第1幅与第960幅二维投影图像的差值图像；

图6是应用本发明发明方法前重建获得的某一三维CT断层图像；

图7是应用本发明发明方法后重建获得的某一三维CT断层图像；

图8是人工设定360度内二维投影图像幅数为900重建获得的某一三维CT断层图像。

具体实施方式

(1)将被扫描物体放置于三维锥束工业CT系统转台，确保任一扫描角度下，物体被锥束覆盖；

(2)以锥束射线对物体实施透照，同时，检台匀速连续旋转，由面阵探测器以固定采样速度连续采集透射过物体的射线投影，获得二维投影图像；

(3)当检台旋转380度时，探测器停止采样，检台和射线源同时停止，即完成一次单圆轨道锥束CT扫描；

(4)当扫描结束，将面阵探测器获取的多个二维投影图像按采集时间顺序排列起来，形成一组二维投影图像I(i，m，n)，其中，i表示扫描角度，m和n是某一探测通道在面阵探测器上的位置；

(5)以如下公式对步骤(4)中二维投影图像进行对数解调，获得一组二维线积分图像：

$P(i,m,n)=\ln \left(\frac{\mathrm{mean}\left(I(i,1:\mathrm{10,1}:K)\right)}{I(i,m,n)}\right)$

(6)以如下公式对步骤(5)中第一幅二维线积分图像与步骤(5)中其他二维线积分图像的相关系数，形成一个一维相关系数数组：

C(i)＝corr(P(1，m，n)，P(i，m，n))

(7)以如下公式搜索步骤(6)获得的一维相关系数数组中的最大值，记录下该最大值在该数组中的序号：

N＝i|_{C(i)＝max(C)}

(8)以如下公式将步骤(4)获得的最大值对应的序号减1，并除2π，结果即为连续扫描三维工业CT的角度增量：

A＝2π/(N-1)

为验证本发明，在实际锥束工业CT系统上进行了实验。实验的具体步骤如下：

(1)设定实验条件。本实验采用的物体是CT空间分辨率测试卡，基体材料为有机玻璃，线对特征为铝；X射线管电压120KV，管电流为1.5mA，射线源到探测器距离1200mm，探测器采样帧频为2帧/秒，转台旋转速度为0.75度/秒，采集的二维投影图像总数1200幅，投影图像尺寸为500*750。

(2)根据设定参数，转台旋转，探测器采集数据，生成一组二维投影图像序列。

(4)根据本发明对数解调公式，对二维投影图像进行解调，生成一组二维线积分图像序列。

(5)根据本发明相关系数计算公式，计算第一幅二维线积分图像与其他所有二维线积分图像的相关系数，形成一维相关系数数组。

(6)搜索一维相关系数数组中的最大值，记录下该最大值在该数组中的序号。

(7)根据本发明角度增量计算公式，计算角度增量。

(8)将角度增量带入FDK三维CT重建算法重建CT图像。

按照实验设定的转台旋转速度和探测器采样帧频，探测器在360°内采集的投影图像的总数为960幅，角度增量为0.375°；按照本发明方法，探测器在360°内采集的投影图像的总数为902幅，角度增量为0.3996°。图4是本发明方法计算得到的一维相关系数数组的曲线图。图中左箭头对应角度为0°时的相关系数，图中右箭头对应角度为360°时的相关系数。这表明探测器在360°内采集的投影图像的总数确实为902幅。图5(a)是第1幅二维投影图像；图5(b)是第902幅二维投影图像；图5(c)是第960幅二维投影图像；图5(d)是第1幅与第902幅二维投影图像的差值图像；图5(e)是第1幅与第960幅二维投影图像的差值图像。与图5(d)比较，图5(e)存在较大的残差。这说明探测器在360°内采集的投影图像的总数确实为902幅，角度增量为0.3996°。图6是以0.375°为角度增量获得的CT重建图像，图7是以0.3996°为角度增量获得的CT重建图像。由图7可知，应用本发明方法后，CT图像伪影和模糊显著降低。为考察结果的准确性，我们以人工设定探测器在360°内采集的投影图像的总数为901幅，角度增量为0.4°，然后进行CT重建，其重建图像如图8所示。比较图7和图8可知，图8出现了较为明显的特征扭曲和模糊。这说明本发明方法确定的角度增量是准确的，同时也从实验上说明角度增量的微小误差确实会在CT图像上导致可见的伪影和模糊。

由图6、7、8可知，本发明方法能正确确定连续扫描三维锥束CT的角度增量，精度高，实现过程简单、高效，不需要特殊的硬件和额外的锥束扫描投影数据，易于工程实现。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种连续扫描三维锥束工业CT角度增量确定方法，其特征在于包括下列步骤：

(1)进行单圆轨道锥束CT扫描，获得一组二维投影图像；

(2)对步骤(1)中所述二维投影图像进行对数解调，获得一组二维线积分图像；

(3)计算步骤(2)中第一幅二维线积分图像与步骤(2)中其他二维线积分图像的相关系数，形成一个一维相关系数数组；

(4)搜索步骤(3)获得的一维相关系数数组中的最大值，记录下该最大值在该数组中的序号；

(5)将步骤(4)获得的最大值对应的序号减1，并除2π，结果即为连续扫描三维工业CT的角度增量。

2.根据权利要求1所述的一种连续扫描三维锥束工业CT角度增量确定方法，其特征在于：所述步骤(1)中进行单圆轨道锥束CT扫描，获得一组二维投影图像步骤为：

(1.1)将被扫描物体放置于三维锥束工业CT系统转台，确保任一扫描角度下，物体被锥束覆盖；

(1.2)以锥束射线对物体实施透照，同时，检台匀速连续旋转，由面阵探测器以固定采样帧频连续采集透射过物体的射线投影，获得二维投影图像；

(1.3)当检台旋转380度时，探测器停止采样，检台和射线源同时停止，即完成一次单圆轨道锥束CT扫描；

(1.4)当扫描结束，将面阵探测器获取的多个二维投影图像按采集时间顺序排列起来，形成一组二维投影图像I(i，m，n)，其中，i表示扫描角度，m和n是某一探测通道在面阵探测器上的位置。

3.根据权利要求1所述的一种连续扫描三维锥束工业CT角度增量确定方法，其特征在于：所述的步骤(2)中对步骤(1)中所述二维投影图像进行对数解调，获得一组二维线积分图像的方法如下：

$P(i,m,n)=\ln \left(\frac{\mathrm{mean}\left(I(i,1:\mathrm{10,1}:K)\right)}{I(i,m,n)}\right)$

其中，P(i，m，n)表示第i个扫描角度对应的二维线积分图像；ln表示自然对数运算；mean表示二维均值运算；1:10表示从1取到10，1:K表示从1取到K，K为二维线积分图像的高度，(1:10，1:K)定义了二维图像上的一个宽度为10高度为K的区域。

4.根据权利要求1所述的一种连续扫描三维锥束工业CT角度增量确定方法，其特征在于：所述的步骤(3)计算步骤(2)中第一幅二维线积分图像与步骤(2)中其他二维线积分图像的相关系数，形成一个一维相关系数数组的方法为：

C(i)＝corr(P(1，m，n)，P(i，m，n))

其中，C(i)表示第1个二维线积分图像与第i个二维线积分图像的相关系数；corr表示二维相关运算；P(1，m，n)表示第1个扫描角度对应的二维线积分图像；P(i，m，n)表示第i个扫描角度对应的二维线积分图像。

5.根据权利要求1所述的一种连续扫描三维锥束工业CT角度增量确定方法，其特征在于：所述的步骤(4)中搜索步骤(3)获得的一维相关系数数组中的最大值，记录下该最大值在该数组中的序号的方法为：

N＝i|_{C(i)＝max(C)}

其中，N表示与第1个二维线积分图像具有最大相关性的二维线积分图像对应的扫描角度序号，也表示360°内探测器采集的图像的个数；C表示第1个二维线积分图像与其他二维线积分图像的相关系数；max表示一维数组取最大值运算。

6.根据权利要求1所述的一种连续扫描三维锥束工业CT角度增量确定方法，其特征在于：所述的将步骤(4)获得的最大值对应的序号减1，并除2π，结果即为连续扫描三维锥束工业CT的角度增量的方法为：

A＝2π/(N-1)

其中，A表示角度增量；N表示与第1个二维线积分图像具有最大相关性的二维线积分图像对应的扫描角度序号。

Images