本发明总体上涉及 CT图像领域， 尤其涉及一种 CT图像重建方法及 装置。 背景技术

自从 1972年 Hounsfield发明了第一台 CT机， CT技术给医学诊断 和工业无损检测带来了革命性的影响， CT已经成为医疗、 生物、 航空 航天、 国防等行业重要的检测手段之一。 随着技术的进步， CT扫描模 式和成像方法也在不断地改进，三维锥束 CT已经成为研究和应用的主 流。 X射线锥束 CT已经在医学临床、 安全检查、 无损检测等领域得到 了广泛的应用。 特别是基于圆轨道扫描的锥束 CT系统由于其在机械、 电控等方面相对简单， 易于工程实现， 因此在临床、 安检和工业无损 检测应用非常广泛。 在圆轨道锥束 CT 重建方法中， 应用最广泛的是 1984年由 Feldkamp等人提出的 FDK方法 ( Feldkamp.L.A., L.C.Davis, and J.W.Kress. Practical cone-beam algorithm [J]. Journal of the Optical Society of America, 1984, (1):612-619. ) 。

FDK方法可以看作是扇束 FBP ( Filtered Backprojection, 滤波反投 影） 方法在三维情况下的近似扩展。 该 FDK方法 ¾括如下步骤： 首先 对投影数据进行加权处理； 然后对不同投影角度的投影数据进行水平 方向上的一维滤波； 最后沿 X射线逆方向进行三维反投影得到物体最 终的三维重建图像。

因此可知 FDK 方法的体素重建值是在 360·投影角度范围内通过 该体素的射线贡献之和，所以圆轨道锥束 FDK方法作为一种近似方法， 具有数学公式简单、 方法计算快速、 易于工程实现等特点， 并且当锥 角比较小的时候 （一般在 ± 5°以内） ， 能够取得很好的重建效果， 因 此在实际工程应用中得到了广泛认可。

但是， FDK方法也存在一定的问题： 由于圆轨道扫描本身不满足 锥束精确重建的数据完备性条件， 存在雷当 （Radon )数据缺失问题； 因此， 当锥束锥角增大时， FDK方法重建结果会产生严重的锥束伪影， 在远离扫描轨道平面内 FDK重建数值快速下降， 使得该方法.在平板探 测器 CT成像系统中的应用受到了很大限制。

为了改善大锥角下圆轨道 FDK方法的图像重建质量， 在 FDK方 法的基础上，提出了多种改进的 FDK方法，例如 P-FDK( parallel FDK )、 T-FDK( tent-FDK )、 HT-FDK( hybrid tent-FDK )、 EFDK( extended FDK ) 等。 在这些改进的 FDK类方法中， 由于 P-FDK和 T-FDK方法简单， 易于在工程中实现， 因此得到了较为广泛的应用， 下面简单介绍一下 这两种方法。

P-FDK 方法是将锥束投影数据通过重排获得平行的扇 束投影数 据， 然后再通过滤波反投影方法重建出物体的三维 图像。 如图 1 ( a ) 所示， 其是使用平板探测器的圆轨道锥束 CT系统扫描示意图， SC ^表 示 X光源在圆轨道上的位置， 表示投影在圆轨道上的角度采样位置。

( ?， _α ，δ)表示平板探测器上的投影数据， （"， W是定义在平板探测器上的 直角坐标系， 用于表示每条 X射线投影在探测器上的位置坐标， R是 圆轨道半径； 图 （b )是利用 P-FDK方法把锥束投影重排成平行束的扇 束投影， 黑色实心点为 X光源， 在中心虚拟探测器上， 重排后的投影 数据由于同一角度下平行扇束到虚拟探测器的 距离不相同， 因此原平 板探测器上一行投影对应到该虚拟探测器不再 是水平的， 而是在一段 沿虚拟探测器中心行外凸的曲线上。 以平板探测器为例， P-FDK 的数 据重排公式如公式（ 1 )。 本发明后续推导过程也皆以平板探测器为例， 其他探测器类型例如柱面探测器等， 都可以根据平板探测器做相应变 形得到， 不做赘述。 ( 1 )

其中， 表示经过 P-FDK重排完后， 其重排投影数据在角度方向 的采样位置。 （^)是经过重排后在中心虚拟探测器上的直角� ��标系， 表 示重排后每条 X射线在虛拟探测器上的位置坐标。

P-FDK方法和 FDK方法的不同之处， 只是由于 P-FDK方法先重 排成平行的扇束使得在反投影时省去了计算加 权系数的过程， 在图像 质量上和 FDK方法没有区别。 不同于 P-FDK, 2000年由 Grass等人提 出的 T-FDK方法在 P-FDK方法基础上做了改进， T-FDK除了将锥束 投影重排成平行的扇束以外， 还在竖直扇束平面内进行二次重排， 即

T-FDK 需要在水平和竖直两个方向上对投影数据进行 重排， 最终使得 T-FDK和 P-FDK的差别在于 T-FDK对投影数据的滤波方向是沿中心 虛拟探测器的水平方向进行的， 如图 1 ( c ) 所示， 而不是沿外凸的曲 线方向。 T-FDK方法的数据重排公式如下，

P ^T ~ ^FDK (θ, t, s) = Ρ(θ - arcsin -, . ^tR , _Ύ ) ( 2 ) 其中， 表示经过 T-FDK重排完后， 其重排投影数据在角度方向 的釆样位置。 （^)是经过 T-FDK重排后在中心虛拟探测器上的直角坐 标系， 表示重排后每条 X射线在虛拟探测器上的位置坐标。

T-FDK相比较 FDK来说， 一方面类似于 P-FDK由于重排成了平 行的扇束省去了反投影时的加权系数， 因此在计算方面更加有效率； 同时， 由于 T-FDK方法对投影数据的滤波是沿中心虚拟探测� ��的水平 方向进行的， 减轻了锥角增大时导致的锥束伪影， 提高了图像重建的 质量， 使得利用大面积平板探测器通过圆轨道扫描实 现大物体的准确 三维成像成为可能。 另外还有一类改进 FDK方法的思路是利用圆轨道 扫描投影中的共轭射线， 通过对共轭投影选取不同的反投影权重系数 来改善重建图像质量， 也取得了较好了效果。

随着平板探测器的日益普及， 使用大面积平板探测器的新型锥束 CT系统越来越多， 对大锥角圆轨道锥束 CT图像重建方法的要求也越 来越高。以目前在牙科疾病诊断中应用较为广 泛的牙科锥束 CT设备为 例， 目前多个厂家的三维牙科 CT使用了 20cm X 25cm平板探测器， X 光源到探测器的距离 设为 70cm , 此时圆轨道扫描对应的锥角大小为 ± 8.13。,由于医生主要依靠图像的 CT数即重建图像像素的数值大小进 行疾病诊断， 因此， 医用 CT对图像的重建数值要求非常高； 如此大锥 角的圆轨道扫描已经远远超出了 FDK方法能够重建的范围， T-FDK重 建结果也没法让人满意。 并且， 随着探测器技术的飞速发展， 更大面 积的平板探测器已经在临床得到了应用， 例如 30cm x 40cm、 43cm χ 43cm 等多种不同尺寸的平板探测器都已经在临床 DR ( Digital Radiography System, 数字式 X射线摄影技术）成像中应用了。 这些更 大面积的探测器能够大大提高锥束 X射线的有效探测面积， 扩大成像 视野，最重要的是能够减少甚至消除锥束投影 数据截断导致的 CT数值 无法准确重建的问题； 因此, 大面积平板探测器在目前和未来的三维

CT影像设备中具有极为广阔的应用空间。 但是， 随着平板探测器面积 的增大， 三维 CT系统的锥角也相应增大， 如何消除大锥角下的严重的 锥束伪影成为一个必须解决的难题。 发明内容

本发明要解决的主要技术问题是提供一种能够 消除大锥角下严重 锥束伪影的 CT图像重建方法及装置。

为了解决上述问题， 本发明 CT 图像重建方法的技术方案包括步 骤：

选取与扫描圆轨道近似曲率的曲线上的同高度 的投影数据； 对所选取的投影数据进行加权处理；

对经过加权处理的投影数据沿水平方向进行滤 波；

对滤波后的投影数据沿射线方向进行三维反投 影。

所述步骤选取与扫描圆轨道近似曲率的曲线上 的同高度的投影数 据包括按下列公式选取投影数据：

P ^C ~ ^FDK (Θ, t, c) = Ρ(θ

其中， ― )表示所选取的投影数据； 表示投影方向； t表示 平行扇束之间的距离； c表示在 Z轴方向上的角度釆样间隔； R表示圆 轨道的半径。

所述步骤对所选取的投影数据进行加权处理包 括按下列公式对所 选取的投影数据进行处理：

其中， φ, ) 表示经过加权处理的投影数据。

所述步骤对经过加权处理的投影数据沿水平方 向进行滤波包括按 下列公式进行滤波： g ^c - ^FDK (0,t,c) = P ^c - ^FDK (e,t,c) h(t)

其中， g ^e — ""^，^)表示经过滤波的投影数据； ®表示卷积； 表 示滤波函数。

所述步骤对滤波后的投影数据沿射线方向进行 三维反投影包括按 下列公式进行三维反投影：

f ^c ~ ^FDK ( , y,z)= \ g ^c ~ ^FDK (Θ, t(x, y, Θ), c(x, y, z, θ)) άθ

0

其中， 皿^，^)表示在 χ轴, y轴 , _z 轴方向上的重建图像。

相应地， 本发明 CT图像重建装置的技术方案包括：

重排单元， 用于选取与扫描圆轨道近似曲率的曲线上的同 高度的 投影数据；

加权单元， 用于对所选取的投影数据进行加权处理；

滤波单元， 用于对经过加权处理的投影数据沿水平方向进 行滤波； 反投影单元， 对滤波后的投影数据沿射线方向进行三维反投 影。 所述重排单元按下列公式选取投影数据：

P ^c ' ^FDK (0,i,c) = P(e

^C-FDK

其中， P - ( 'c)表示所选取的投影数据； 表示投影方向； 表示 平行扇束之间的距离； c表示在 ζ轴方向上的角度采样间隔； R表示圆 轨道的半径。

所述加权单元按下列公式对所选取的投影数据 进行处理：

其中， ? ^c -舰 (θ, )表示经过加权处理的投影数据 ,

所述滤波单元按下列公式进行一维斜坡滤波： g"-^ (θ, t, c) = P ^C - ^FDK φ, t, c) ® hit) ： ^ p^-^^e _{y tj C} . h t - t ^r )dr 其中， 皿 ( ， ^e )表示经过滤波的投影数据； ®表示卷积； 表 示滤波函数。

所述反投影单元按下列公式进行三维反投影：

f- ^FDK {x, y, z) = ζ, θ))άθ

其中， 皿 <， ^ )表示在 X轴, y轴 , _Z 轴方向上的重建图像。

与现有技术相比， 本发明 CT图像重建方法及装置的有益效果为： 首先， 本发明对投影数据进行重排是采用选取与扫描 圆轨道近似 曲率的曲线上的同高度的投影数据， 从而在平行于 Z轴的扇束平面内 沿虚拟中心探测器向中心行内凹的曲线上采样 ， 这样大大提高了在大 锥角下重建方法的数值精度， 有效地抑制了大锥角带来的锥束伪影。

另外， 本发明执行效率高， 稳定性强。 附图说明

为了对本公开内容更透彻的理解， 下面参考结合附图所进行的下 列描述， 在附图中：

图 1(a)是平板探测器圆轨道锥束 CT扫描示意图；

图 1(b)是采用 P-FDK方法对图 1(a)所得投影数据重排后的示意图； 图 1(c)是采用 T-FDK方法对图 1(a)所得投影数据重排后的示意图； 图 1(d)是采用本发明 CT图像重建方法对图 1(a)所得投影数据重排 后的示意图；

图 2是图 l(a)、 1(b), 1(c) 和 1(d)沿旋转轴 Z轴的俯视图； 图 3是图 2过 SC的竖直切面的侧视图；

图 4是本发明 CT图像重建方法的流程图；

图 5(a)示出了圆轨道 CT扫描所采集的投影数据， 其是由在各个角 度下平板探测器中心层数据组成的正弦图；

图 5(b)是经过本发明 C-FDK方法对投影数据进行重排后， 各个角 度下虚拟探测器中心层数据组成的正弦图； 图 5(c)是对图 5(b)所示的投影数据进行加权处理后所得到的� �影 数据的示意图；

图 6(a) 是三维 Shepp-logan头模型在一竖直截面的准确图像的示 意图；

图 6(b)是釆用 FDK方法所得到的重建结果的示意图 ·；

图 6(c)是采用 T-FDK方法所得到的重建結果的示意图；

图 6(d)是采用本发明 C-FDK方法所得到的重建结果的示意图； 图 7是选取沿图 6(a)、 图 6(b) 、 图 6(c) 和图 6(d)截面的竖直中心 的剖面线后重建结果的示意图；

图 8是本发明 CT图像重建装置的示意图。 具体实施方式

下面将详细描述本发明的具体实施例， 但本发明并不限于下述具 体实施例。

如图 4所示， 本发明 CT图像重建方法包括步骤：

1 )选取与扫描圆轨道近似曲率的曲线上的同高� �的投影数据；

2 )对所选取的投影数据进行加权处理；

3 )对经过加权处理的投影数据沿水平方向进行� �波；

4 )对滤波后的投影数据沿射线方向进行三维反� �影。

从上述可知，本发明 CT图像重建方法进行重排的方式是选取与扫 描圆轨道近似曲率的曲线上的同高度的投影数 据， 接着对所选取的投 影数据进行加权处理， 然后将经过加权处理的投影数据沿水平方向进 行滤波， 最后再对滤波后的投影数据沿射线方向进行三 维反投影即可 得到被扫描物体的三维 CT图像。

对于与扫描轨道曲率近似的曲率而言，假设扫 描圆轨道曲率为 1/R , 该近似平均曲率范围可以是 1/ 2R〜2R。对于滤波，可以采用本领域技术 人员可以获知的任何滤波方法， 例如， 滤波核可以是最基本的斜坡滤 波核， 也可以是对标准斜坡滤波在频域进行过平滑处 理后的滤波核， 比如常用的 S-L滤波核等 ( S-L滤波核由 L.A.Shepp和 B.F.Logan于 1974 年提出） 。

本说明书中是以平板探测器为例来描述本发明 的技术方案的。 当 然， 本发明可以应用到诸如柱面探测器之类的其它 面阵探测器。 如图 1 ) 所示， 首先， 定义平面内锥束 X射线源的扫描路径为 一个圆周 C^ - R^osAsin ) ; R是圆轨道半径； 表示光源点对应的角 度参量； O是坐标原点和圆轨道的中心， 也就是旋转中心。 ( ?，α, έ)表 示 X射线照射被扫描物体后在平板探测器 1上采集的投影数据， 其中 (α, >)表示二维平板探测器 1上某个投影点位置的横纵坐标。

根据 P-FDK方法的性质， 该方法对锥束投影数据进行重排后选取 的是在中心虛拟探测器 2上在外凸的曲线上的投影数据， 如图 1 ( b ) 所示。 下面请参见图 2 , 图 2是沿旋转轴 Z轴从上往下的俯视图， 此时 面阵列探测器成为一条直线或曲线， 圆轨道上的实心点 S表示 X射线 光源。 在图 2中， P-FDK方法重排数据的过程可以看成是将图 2中经 过 OP的曲线上的同高度投影数据进行组合， 即得到 P-FDK的数据， 也 就是公式（ 1 )的结果， 其中 OP曲线即为由 P-FDK重排公式（ 1 )定义 的一条曲线。 而 T-FDK方法重排数据的过程可以看成是将图 2中经过 OO'的直线上的同高度投影数据进行组合， 即得到了 T-FDK的数据， 也就是公式（2 ) 的结果， 其中 OO'曲线即为由 T-FDK重排公式 （2 ) 定义的一条曲线。

而本发明 CT 图像重建方法 （简称 C-FDK ) 进行数据重排的过程 则可以看成是将图 2中经过 OC的曲线上的同高度投影数据进行组合， 即得到了 C-FDK的数据， 本例中 OC曲线是以半径为 R的一段圆弧， 也 就是和扫描圆轨道相同曲率的一段弧线， 且 OC的曲率中心在过 O点且 平行于 SC的直线上。 当然， 0C曲线的曲率可以不准确等于扫描圆轨道 的曲率， 可以在圆轨道曲率的 50%〜200%范围之内， 也即在此范围内 的投影数据同样可以达到近似本发明的效果。

本发明 C-FDK方法对投影数据的重排是在中心虚拟探测� ��的水平 和竖直两个方向上完成， 在水平方向与 P-FDK以及 T-FDK是相同的； 在竖直方向上， 如图 3所示， 对重排后的数据： , , 、 c = b K 5 )

其中， c表示 C-FDK方法重排后中心虚拟探测器沿竖直方向的� �� 标， S表示角度为 且到中心射线距离为 t的 X射线与轨道的交点位置， 因此， SC表示位于 S处的 X射线源点沿该条 X射线到曲线 OC的距离。 SP 表示位于 S处的 X射线源点沿该条 X射线到曲线 OP的距离。

图 3是沿图 2中过 SC的竖直切面， 其中， 与 SC所在的轴成 角的 线为 X射线。 根据图 2中 OC弧线的定义， 图 3 中 S (：、 SP的长度都可 以通过空间几何关系计算得到， 其具体长度如下：

因此， 本发明 C-FDK方法中对圆轨道锥束 CT投影数据的重排公 式:^下： ( 7 )

其中， 表示投影方向， t表示平行扇束之间的距离， c表示在 Z 轴方向上的角度采样间隔。

利用上面公式 （7 ) 完成 C-FDK 方法的数据重排后， 下一步是对 投影数据进行加权处理， 加权系数可以为 cos , 其中 是每条 X射线与 其沿 Z轴方向投影到中心平面上的投影线之间的夹� �， 如图 3所示， 根据空间几何关系：

即对投影进行加权处理的公式为 ( 9 )

然后， 对完成上述处理的投影数据沿水平方向进行一 维斜坡滤波: g ^c - ^mK φ, t, c) = P ^C - ^FDK (θ, t, c) ® h(t)

¾ ( in ) 上式中， ®表示卷积， 而; t)是滤波函数， 一般采用 Ramp滤波器。 最后， 将滤波后的投影数据沿射线方向进行三维反投 影即可重建 出被扫描物体的三维 CT图像， 反投影公式如下：

2π

f ^c — ^FDK (x, y,z)= \ g ^c _ ^FDK (Θ, t(x, y, Θ), c(x, y, z, θ)) άθ ( 11)

0

其中， 探测器上的投影位置计算公式如下 t(x, }>,θ) = y cos ^ - x sin ^ ( 12)

z-{2^R ² -f -R)

c(x,y,z,0) ( 13) 如图 5 (a) 所示， 示出了圆轨道 CT扫描所采集的投影数据， 其 疋由在各个角度下平板探测器中心层数据组成 的正弦图； 图 5 (b) 则 是经过本发明 C-FDK方法对投影数据进行重排后， 各个角度下虚拟探 测器中心层数据组成的正弦图； 图 5 (c) 则是对图 5 (b) 所示的投影 数据进行加权处理后所得到的投影数据。 接着对经过加权处理的投影 数据进行一维斜坡滤波， 最后， 再进行反投影即可得到被扫描物体的 CT图像， 如图 6 (d) 所示。

相应地， 如图 8所示， 本发明 CT图像重建装置包括：

重排单元 1,用于选取与扫描圆轨道近似曲率的曲线上的� ��高度的 投影数据；

加权单元 2， 用于对所选取的投影数据进行加权处理；

滤波单元 3, 用于对经过加权处理的投影数据沿水平方向进 行滤 波；

反投影单元 4, 对滤波后的投影数据沿射线方向进行三维反投 影。 优选地， 所述重排单元 1按下列公式选取投影数据： 其中， P - ( ，c)表示所选取的投影数据； 表示投影方向； t表示 平行扇束之间的距离； c表示在 Z轴方向上的角度釆样间隔； R表示圆 轨道的半径。

优选地， 所述加权单元 2按如下公式对所选取的投影数据进行处 理：

PC— (θ, ί, C) = - ' - . pC-FDK _φ , _t , _c)

V5R ² - 4t ² - 4R /R ² - t ² + c ²

其中， ( ， c)表示经过加权处理的投影数据。

优选地， 所述滤波单元 3按如下公式进行滤波： gd (Θ, t, c) = P ^C - ^FDK (Θ, t, c) ® hit)

其中， _ ( ， 表示经过滤波的投影数据； ®表示卷积； 表 示滤波函数。

优选地， 所述反投影单元 4按如下公式进行三维反投影：

其中， /—皿(^ 表示在 X轴, y轴 ,ζ轴方向上的重建图像。

由于本发明 CT 图像重建装置的技术方案所包含的技术特征与 本 发明 CT重建方法的技术方案所包含的技术特征相对� ��，因此不再对本 发明 CT图像重建装置的技术方案进行详细描述。

本申请的发明人利用三维 Shepp-Logan 头模型进行了数值模拟实 验，并与 FDK方法和 T-FDK方法进行了实验对比来验证本发明的技术 方案。 在数值模拟实验中， 三维头模型被限制在 lmm的球体内， 模型 中心即为 CT扫描的旋转中心，. X光源到旋转中心距离为 4mm, X光 源到探测器的距离为 8mm, 平板探测器尺寸为 4mmx4mm, 探测单元 数量为 256x256 , 在 360度范围内角度均匀釆集 360个锥束投影， 然后 进行三维 CT图像重建。 根据上述几何定义， 可以计算出在此实验中， 锥束 X射线的最大锥角约为 14度， 即本次实验的 X射线锥角范围为 ±14°。使用上述扫描条件下的圆轨道锥束 CT投影数据，分别采用 FDK 方法、 T-FDK方法和本发明 C-FDK方法进行三维 CT图像重建。

图 6 ( a )是三维 Shepp-logan头模型在一竖直截面的准确图像， 所 述竖直截面是三维头模型在 y=-0.25处，沿平行于 x-z平面予以选取的； 图 6 ( b )是采用 FDK方法所得到的重建结果； 图 6 ( c )是采用 T-FDK 方法所得到的重建结果； 图 6 ( d )是采用本发明 C-FDK方法所得到的 重建结果。

从上述所得到的各重建图像可以看出： 本发明 C-FDK方法能够在 ±14。条件下， 很好地重建出三维头模型的图像， 克服了现有方法在大 锥角圆轨道锥束 CT重建中普遍存在的锥束伪影难题，很好地解� ��了大 锥角圆轨道 CT图像重建问题。

为了更准确地分析上述重建结果的数值准确性 ， 图 7 选取沿图 6 截面的竖直中心的剖面线， 其中， 横坐标表示该剖面线上点的长度坐 标， 纵坐标表示该剖面线上点的线性衰减系数， 长虚线表示 FDK方法 的重建结果， 点划线表示 T-FDK方法的重建结果， 短虚线表示本发明 C-FDK方法的重建结果， 实线表示头模型的准确数值。 从该剖面线可 以更清楚地看到：现有的 FDK和 T-FDK方法当锥角增大时其重建数值 呈快速下降趋向， 离模型的准确值越来越远， 而本发明 C-FDK方法即 使在 ± 14°锥角范围内仍然能够较为准确地重建出原� �模型的准确值， 基本解决了大锥角圆轨道 CT图像重建的难题。

最后， 需要特别说明的是： 本发明 C-FDK方法的核心之处在于沿 中线虛拟探测器内凹的方向上进行投影数据的 滤波， 在上述推导过程 中， 本发明是将 X射线和图 2中 OC曲线上同高度的投影重排为同一高 度 c的 C-FDK数据 P ^c -皿 (6>, t, _c )。 而 OC曲线则是和圆扫描轨道弧度相同 的一段圆弧， 此处 OC曲线的选取不是唯一的， 可以在该曲率附近做小 范围的改变， 也同样应该在本发明专利的权利要求范围内。

虽然上述已经结合附图描述了本发明的具体实 施例， 但是本领域 技术人员在不脱离本发明的精神和范围的情况 下， 可以对本发明进行 各种改变、 修改和等效替代。 这些改变、 修改和等效替代都意为落入 随附的权利要求所限定的精神和范围之内。