CN101825433A - 扇束2d-ct扫描***转台旋转中心偏移量的测量方法 - Google Patents

Abstract

本发明公开了一种扇束2D-CT扫描***转台旋转中心偏移量的测量方法，是将扇束投影重排为平行束投影，根据投影角相差π的平行束投影相似的原理，将重排后的平行束正弦图中[0°180°]和[180°360°]对应的投影信息取出，获得第一正弦图和第二正弦图；将第二正弦图水平翻转后形成的第三正弦图与第一正弦图进行互相关运算，从而求出转台旋转中心偏移量τ₀。本发明的方法不需要制作专门的校正试件，直接利用被扫描断层的正弦图进行测量。有效保证了扇束2D-CT***转台旋转中心偏移量的测量精度，提高了重建图像的质量。

Description

扇束2D-CT扫描***转台旋转中心偏移量的测量方法

技术领域

本发明涉及一种适用于扇束2D-CT(Two Dimensional ComputedTomography)扫描***转台旋转中心偏移量的测量，可用于医学和工业领域射线数字成像(DR-Digital Radiography)、二维、三维计算机断层扫描(2D-CT/3D-CT)成像过程中的相关测量。

背景技术

对于基于线阵列探测器的扇形束射线扫描2D-CT，其扫描原理如图1所示。转台4位于射线源焦点1与线阵列探测器6之间；转台4处于扇束射线3的辐照场中；中心射线2垂直于线阵列探测器6，其垂足记为S₀，S₀也是线阵列探测器6的中心位置的坐标；根据2D-CT重建算法的要求，中心射线2要通过转台4的旋转中心点O并垂直于线阵列探测器6。然而，对于实际的扇束2D-CT***，由于机械安装上的误差，会使得旋转中心点O偏移中心射线2，这种偏移量τ₀的存在会影响CT图像的重建质量，导致重建伪影的出现。为了得到精确的重建结果，必须准确测量出旋转中心点O相对于中心射线2的偏移量τ₀，然后将获取的偏移量τ₀作为转台4的中心校正参数输入图像重建算法中。

发明内容

本发明的目的是提出一种适用于扇束2D-CT扫描***转台旋转中心偏移量的测量方法，是将扇束投影重排为平行束投影，根据投影角相差π的平行束投影相似的原理，将重排后的平行束正弦图中[0°180°]和对应的投影信息取出，获得第一正弦图和第二正弦图；将第二正弦图水平翻转后形成的第三正弦图与第一正弦图进行互相关运算，从而求出转台旋转中心偏移量τ₀。

本发明是一种适用于扇束2D-CT扫描***转台旋转中心偏移量的测量方法，该测量方法包括有下列实施步骤：

步骤一：在2D-CT扫描***中，线阵列探测器6和射线源焦点1不动，让转台4绕旋转中心点O作圆周转动；

在本发明中，转台4绕旋转中心点O的转动，使得线阵列探测器6能够采集到不同投影角下的被扫描断层5的投影数据信息，将所有投影角下的投影数据合成为一幅扇束射线扫描正弦图p^fan(θ，s)。关于转动角度(也是投影角)的设置可以精确到每转0.1°～1°进行一次投影数据的采集。

步骤二：将扇束射线扫描正弦图p^fan(θ，s)重排为平行束扫描正弦图p(θ，s)，p(θ，s)为二维矩阵，该矩阵的每一行代表每一投影角下的投影数据，θ＝β时的投影数据记为p_β(s)，p_β(s)为p(θ，s)的某一行；

步骤三：将正弦图p(θ，s)中θ∈[0°180°]范围内的对应投影数据取出合成第一正弦图p_β ¹(s)；将正弦图p(θ，s)中θ∈[180°360°]范围内的对应投影数据取出合成第二正弦图p_β ²(s)；可以看出，p_β ¹(s)和p_β ²(s)均为二维矩阵，其高度为正弦图p(θ，s)的1/2，p_β ¹(s)为p(θ，s)的前半部分，p_β ²(s)为p(θ，s)的后半部分；

步骤四：将p_β ²(s)进行水平翻转，得到第三正弦图p_β ³(s)，p_β ³(s)与p_β ²(s)的关系为：

N为线阵列探测器探测单元的总数目。

步骤五：将第一正弦图p_β ¹(s)与第三正弦图p_β ³(s)的对应行进行互相关运算，那么对应行的互相关函数最大值对应的坐标就是p_β ¹(s)与p_β ³(s)对应行之间的位移差数列记为d(i)，i＝1～M，M表示p_β ¹(s)或p_β ³(s)的高度。

步骤六：求数列d(i)的均值τ′，该均值τ′的1/2等于转台旋转中心点的偏移量τ₀，即 ${\mathrm{\τ}}_0=\frac{1}{2M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}d\left(i\right).$

采用本发明的测量方法进行扇束2D-CT扫描***转台旋转中心偏移量测量，具有如下的优点：(1)不需要制作专门的校正试件，直接利用被扫描断层的正弦图进行测量。有效保证了扇束2D-CT***转台旋转中心偏移量的测量精度，提高了重建图像的质量。(2)将互相关方法应用到扫描正弦图的解算过程中，提高了本发明测量方法的抗噪能力。(3)测量算法简单，可作为扇束2D-CT***的辅助校验模块，自动监测转台旋转中心的变化。

附图说明

图1是基于线阵列探测器的扇形束射线2D-CT扫描原理图。

图1A是将图1中的扇形束投影重排为平行束投影的扫描原理图。

图2是θ＝β时平行束射线扫描原理图。

图2A是θ＝β+π时平行束射线扫描原图。

图3是p_β(s)、p_β+π(s)与p′_β+π(s)灰度曲线图。

图4是重排为平行束射线的正弦图。

图4A是第一正弦图。

图4B是第二正弦图。

图4C是第三正弦图。

图4D是第三正弦图与第一正弦图的互相关图像。

图5A与图5B是中心偏移量τ₀校正前后的重建图像。

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

参见图1所示，扇束2D-CT***中的转台4位于射线源焦点1与线阵列探测器6之间；转台4处于扇束射线3的辐照场中；中心射线2垂直于线阵列探测器6，其垂足记为S₀，S₀也是线阵列探测器6的中心位置的坐标。转台4的坐标系为xoy。

线阵列探测器6和射线源焦点1不动，转台4绕旋转中心点O在360°范围内旋转(即转台4绕旋转中心点O作圆周运动)，被扫描断层5在转台4的拖动下同步旋转，扇束射线3穿过被扫描断层5后强度会发生变化，线阵列探测器6在转台4每转过一个固定角度(一般为0.1°～1°)的位置采集射线强度信息。

线阵列探测器6由多个独立的探测单元排列成一行，射线强度信息被每个探测单元采集并经过A/D转换传输到计算机(计算机是一种能够按照事先存储的程序，自动、高速地进行大量数值计算和各种信息处理的现代化智能电子设备。最低配置为CPU 2GHz，内存2GB，硬盘180GB；操作***为windows2000/2003/XP。)，每个探测器单元采集的数据称为投影值。线阵列探测器6每次采集到一行数据，将其在360°范围内采集到的所有行数据依次排列，便构成二维矩阵，称为正弦图。扇束射线扫描正弦图记为p^fan(θ，S)，θ表示转台的旋转角度，也代表投影角，S代表探测器单元坐标。

根据庄天戈公开的《CT原理与算法》可知，扇束射线投影数据可以重排为平行束射线投影数据，如图1A所示。本发明的测量方法是基于将扇束投影数据重排成的平行束投影数据进行转台旋转中心偏移量τ₀的精确测量。在文献(庄天戈，CT原理与算法，上海交通大学出版社，1992)中，公开了基于扇束投影数据重排为平行束投影数据的图像重建算法，其基本思想为：根据扇形射线束扫描的几何关系，利用数学插值方法，将扇束投影数据重排为平行束投影数据，再对重排后的平行束投影数据进行滤波，利用平行束的投影几何关系计算投影地址，最后进行反投影运算得到重建图像。

设转台4的坐标系为xoy，被扫描断层5的中心点O′(也称被扫描断层中心点O′)与旋转中心点O的连线O′O与x轴的夹角记为投影角θ，且θ∈[0°360°]，θ也代表转台相对于x的转角。

当投影角θ的角度取值为β，β∈[0°180°](θ＝β)时，线阵列探测器6采集到的被扫描断层5的投影数据记为p_β(s)，或称第一投影信息p_β(s)，该数据的灰度曲线如图3(a)中曲线所示。

当投影角θ的角度取值在π+β时(θ＝β+π)，线阵列探测器6采集到的被扫描断层5的投影数据记为p_β+π(s)，或称第二投影信息p_β+π(s)，该数据的灰度曲线如图3(b)中曲线所示。

对比不同投影角度(θ＝β或者θ＝β+π)时的投影灰度曲线图，可知第一投影信息p_β(s)和第二投影信息p_β+π(s)存在着内在的联系为：

图3(a)表示第一投影信息p_β(s)的灰度曲线，坐标轴O_dS表示线阵列探测器6的长度方向的坐标轴，线阵列探测器6的长度记为N，这也代表了线阵列探测器探测单元的总数。Q₁为灰度曲线p_β(s)的最高点，其在坐标轴O_dS上对应坐标为S₂，线阵列探测器6的中心位置在坐标轴O_dS上对应坐标为S₀，S₀也是中心射线2垂直于线阵列探测器6的垂足，转台4的旋转中心点O在坐标轴O_dS上对应坐标值记为R₀，由图可知，

S₀-R₀＝τ₀，τ₀表示中心射线2偏移转台4的旋转中心点O的偏移量。

图3(b)表示第二投影信息p_β+π(s)的灰度曲线，Q₂为灰度曲线p_β+π(s)的最高点，其在坐标轴O_dS上对应坐标为S₁。

由平行束扫描的投影原理可知，图3(a)和图3(b)的灰度曲线在坐标轴O_dS上以R₀点为中心呈左右对称关系，即Q₁点与Q₂点关于R₀点对称，满足的关系为：

$|S_1-R_0|=|S_2-R_0|=\frac{1}{2}|S_2-S_1|$

即S₁＝S₂-(S₂-S₁)＝S₂-2(S₂-R₀)＝2R₀-S₂

将第二投影信息p_β+π(s)在坐标轴O_dS上进行左右翻转，便得到翻转投影信息p′_β+π(s)，p′_β+π(s)与p_β+π(s)的关系为：

p′_β+π(s)＝p_β+π(N-s)   s∈[0  N]      (2)

图3(c)为翻转投影信息p′_β+π(s)的灰度曲线，Q₃为p′_β+π(s)的最高点，其在坐标轴O_dS上对应坐标为S₃，联立式(1)、式(2)得p′_β+π(s)的投影坐标关系：

S₃＝N-S₁＝N-(2R₀-S₂)    (3)

对比图3(a)与图3(c)可知，p_β(s)与p′_β+π(s)的灰度曲线形状相同，只是在坐标轴O_dS方向上有一定的位移差，该位移差即是点Q₁与Q₃的坐标差记为τ′，其表达式为：

τ′＝S₃-S₂       (4)

将式(3)代入式(4)得：

${\mathrm{\τ}}^{\′}=S_3-S_2=N-(2R_0-S_2)-S_2=N-2R_0$

$=2\×(\frac{N}{2}-R_0)=2\×(S_0-R_0)=2{\mathrm{\τ}}_0$

即： ${\mathrm{\τ}}_0=\frac{{\mathrm{\τ}}^{\′}}{2}.$

由式(4)可以得出结论：利用两个投影角相差180°的投影数据可以求出转台旋转中心偏移量τ₀，其方法为将其中一个投影在坐标轴O_dS上进行左右翻转，求取另一个投影数据与翻转后的投影数据在坐标轴O_dS上的位移差τ′，那么转台旋转中心的偏移量τ₀等于τ′的一半。

本发明是一种适用于扇束2D-CT扫描***转台旋转中心偏移量的测量方法，该测量方法包括有下列实施步骤：

步骤一：在2D-CT扫描***中，线阵列探测器6和射线源焦点1不动，让转台4绕旋转中心点O作圆周转动；

在本发明中，转台4绕旋转中心点O的转动，使得线阵列探测器6能够采集到不同投影角下的被扫描断层5的投影数据信息，将所有投影角下的投影数据合成为一幅扇束射线扫描正弦图p^fan(θ，s)。关于转动角度(也是投影角)的设置可以精确到每转0.1°～1°进行一次投影数据的采集。

步骤二：将扇束射线扫描正弦图p^fan(θ，s)重排为平行束扫描正弦图p(θ，s)，p(θ，s)为二维矩阵，该矩阵的每一行代表每一投影角下的投影数据，θ＝β时的投影数据记为p_β(s)，p_β(s)为p(θ，s)的某一行；

步骤三：将正弦图p(θ，s)中θ∈[0°180°]范围内的对应投影数据取出合成第一正弦图p_β ¹(s)；将正弦图p(θ，s)中θ∈[180°360°]范围内的对应投影数据取出合成第二正弦图p_β ²(s)；可以看出，p_β ¹(s)和p_β ²(s)均为二维矩阵，其高度为正弦图p(θ，s)的1/2，p_β ¹(s)为p(θ，s)的前半部分，p_β ²(s)为p(θ，s)的后半部分；

步骤四：将p_β ²(s)进行水平翻转，得到第三正弦图p_β ³(s)，p_β ³(s)与p_β ²(s)的关系为：

N为线阵列探测器探测单元的总数目。

步骤五：将第一正弦图p_β ¹(s)与第三正弦图p_β ³(s)的对应行进行互相关运算，那么对应行的互相关函数最大值对应的坐标就是p_β ¹(s)与p_β ³(s)对应行之间的位移差数列记为d(i)，i＝1～M，M表示p_β ¹(s)或p_β ³(s)的高度。

步骤六：求数列d(i)的均值τ′，该均值τ′的1/2等于转台旋转中心点的偏移量τ₀，即 ${\mathrm{\τ}}_0=\frac{1}{2M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}d\left(i\right).$

下面对本发明的测量方法进行实验验证：

成像条件：射线源焦点1的尺寸为2.5mm，管电压420kV，管电流2.5mA。线阵列探测器6的探元尺寸为0.083mm，探测单元数目为2048，焦距为2600mm，投影角数目为1800。

在此成像条件下对涡轮叶片进行扫描，得到扇束正弦图，如图4所示。将图4所示的扇束正弦图中的前1/2和后1/2对应的投影信息取出得到第一正弦图(如图4A)和第二正弦图(如图4B)，记为p_β ¹(s)和p_β ²(s)。将p_β ²(s)进行水平翻转得到第三正弦图p_β ³(s)(如图4C)。将第一正弦图p_β ¹(s)与第三正弦图p_β ³(s)的对应行进行互相关运算得到互相关图像(如图4D)，互相关图像的每行灰度最大值对应的坐标记为d(i)，求数列d(i)的均值记为τ′，该值的1/2等于转台旋转中心的偏移量τ₀， ${\mathrm{\τ}}_0=\frac{1}{2M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}d\left(i\right).$

实验采用的扇束2D-CT***原始标定的转台旋转中心偏移量为0，但是使用三年后，转台的旋转中心位置有所变动，此时如果利用原始标定值会造成重建图像出现伪影，图5A为利用原始值的重建结果；利用本发明的测量方法测得该***转台旋转中心的偏移量为-2.9(单位：像素)，利用该值对重建参数进行修正，得到的重建图像如图5B所示。可见，利用本发明的测量方法可以有效提高扇束2D-CT***的图像重建精度，提高CT图像的细节分辨力。

Claims (3)

1.一种扇束2D-CT扫描***转台旋转中心偏移量的测量方法，其特征在于包括有下列测量步骤：

步骤一：在2D-CT扫描***中，线阵列探测器(6)和射线源焦点(1)不动，让转台(4)绕旋转中心点O作圆周转动；所述转台(4)绕旋转中心点O的转动，使得线阵列探测器(6)能够采集到不同投影角下的被扫描断层(5)的投影数据信息，将所有投影角下的投影数据合成为一幅扇束射线扫描正弦图p^fan(θ，s)；

步骤二：将扇束射线扫描正弦图p^fan(θ，s)重排为平行束扫描正弦图p(θ，s)，p(θ，s)为二维矩阵，该矩阵的每一行代表每一投影角下的投影数据，θ＝β时的投影数据记为p_β(s)，p_β(s)为p(θ，s)的某一行；

步骤三：将正弦图p(θ，s)中θ∈[0°180°]范围内的对应投影数据取出合成第一正弦图p_β ¹(s)；将正弦图p(θ，s)中θ∈[180°360°]范围内的对应投影数据取出合成第二正弦图p_β ²(s)；p_β ¹(s)和p_β ²(s)均为二维矩阵，p_β ¹(s)为p(θ，s)的前半部分，p_β ²(s)为p(θ，s)的后半部分；

步骤四：将p_β ²(s)进行水平翻转，得到第三正弦图p_β ³(s)，p_β ³(s)与p_β ²(s)的关系为p_β ³(s)＝p_β ²(N-s)，N为线阵列探测器(6)中探测单元的总数目；

步骤五：将第一正弦图p_β ¹(s)与第三正弦图p_β ³(s)的对应行进行互相关运算，那么对应行的互相关函数最大值对应的坐标就是p_β ¹(s)与p_β ³(s)对应行之间的位移差数列记为d(i)，i＝1～M，M表示p_β ¹(s)或p_β ³(s)的高度；

步骤六：求位移差数列d(i)的均值τ′，该均值τ′的1/2等于转台旋转中心点的偏移量τ₀，即

2.根据权利要求1所述的扇束2D-CT扫描***转台旋转中心偏移量的测量方法，其特征在于：转台(4)绕旋转中心点O每次转动的角度设置为0.1°～1°。

3.根据权利要求1所述的扇束2D-CT扫描***转台旋转中心偏移量的测量方法，其特征在于：对比不同投影角度θ＝β或者θ＝β+π时的投影灰度曲线图，可知第一投影信息p_β(s)和第二投影信息p_β+π(s)存在着内在的联系为：

所述β表示的是投影角θ的取值，β∈[0°180°]；

在第一投影信息p_β(s)的灰度曲线中，坐标轴O_dS表示线阵列探测器(6)的长度方向的坐标轴，线阵列探测器(6)的长度记为N，这也代表了线阵列探测器探测单元的总数；Q₁为灰度曲线p_β(s)的最高点，其在坐标轴O_dS上对应坐标为S₂，线阵列探测器(6)的中心位置在坐标轴O_dS上对应坐标为S₀，S₀也是中心射线(2)垂直于线阵列探测器(6)的垂足，转台(4)的旋转中心点O在坐标轴O_dS上对应坐标值记为R₀，则有

S₀-R₀＝τ₀，τ₀表示中心射线(2)偏移转台(4)的旋转中心点O的偏移量；

在第二投影信息p_β+π(s)的灰度曲线中，Q₂为灰度曲线p_β+π(s)的最高点，其在坐标轴O_dS上对应坐标为S₁；由平行束扫描的投影原理可知，第一投影信息p_β(s)的灰度曲线与和第二投影信息p_β+π(s)的灰度曲线在坐标轴O_dS上以R₀点为中心呈左右对称关系，即Q₁点与Q₂点关于R₀点对称，满足的关系为：

$|S_1-R_0|=|S_2-R_0|=\frac{1}{2}|S_2-S_1|---\left(1\right)$

即S₁＝S₂-(S₂-S₁)＝S₂-2(S₂-R₀)＝2R₀-S₂

将第二投影信息p_β+π(s)在坐标轴O_dS上进行左右翻转，便得到翻转投影信息p′_β+π(s)，p′_β+π(s)与p_β+π(s)的关系为：

p′_β+π(s)＝p_β+π(N-s)   s∈[0 N]           (2)

在翻转投影信息p′_β+π(s)的灰度曲线，Q₃为p′_β+π(s)的最高点，其在坐标轴O_dS上对应坐标为S₃，联立式(1)、式(2)得p′_β+π(s)的投影坐标关系：

S₃＝N-S₁＝N-(2R₀-S₂)                          (3)

对比第一投影信息p_β(s)的灰度曲线与翻转投影信息p′_β+π(s)的灰度曲线，则有p_β(s)与p′_β+π(s)的灰度曲线形状相同，只是在坐标轴O_dS方向上有一定的位移差，该位移差即是点Q₁与Q₃的坐标差记为τ′，其表达式为：

τ′＝S₃-S₂                                  (4)

将式(3)代入式(4)得：

${\mathrm{\τ}}^{\′}=S_3-S_2=N-(2R_0-S_2)-S_2=N-2R_0$

$=2\×(\frac{N}{2}-R_0)=2\×(S_0-R_0)=2{\mathrm{\τ}}_0$

即： ${\mathrm{\τ}}_0=\frac{{\mathrm{\τ}}^{\′}}{2}.$

Images