CN102456227B - Ct图像重建方法及装置 - Google Patents

Abstract

本发明公开了一种CT图像重建方法及装置。本发明CT图像重建方法包括步骤：选取与扫描圆轨道近似曲率的曲线上的同高度的投影数据；对所选取的投影数据进行加权处理；对经过加权处理的投影数据沿水平方向进行滤波；对滤波后的投影数据沿射线方向进行三维反投影。本发明CT图像重建装置包括重排单元、加权单元、滤波单元和反投影单元。采用本发明的技术方案能够有效地消除大锥角下的锥束伪影。

Description

CT图像重建方法及装置

技术领域

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

背景技术

自从1972年Hounsfield发明了第一台CT机，CT技术给医学诊断和工业无损检测带来了革命性的影响，CT已经成为医疗、生物、航空航天、国防等行业重要的检测手段之一。随着技术的进步，CT扫描模式和成像方法也在不断地改进，三维锥束CT已经成为研究和应用的主流。X射线锥束CT已经在医学临床、安全检查、无损检测等领域得到了广泛的应用。特别是基于圆轨道扫描的锥束CT***由于其在机械、电控等方面相对简单，易于工程实现，因此在临床、安检和工业无损检测应用非常广泛。在圆轨道锥束CT重建方法中，应用最广泛的是1984年由Feldkamp等人提出的FDK方法(Feldkamp.L.A.，L.C.Davis，andJ.W.Kress.Practical cone-beam algorithm[J].Journal of the OpticalSociety 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***扫描示意图，S(β)表示X光源在圆轨道上的位置，β表示投影在圆轨道上的角度采样位置。P(β，a，b)表示平板探测器上的投影数据，(a，b)是定义在平板探测器上的直角坐标系，用于表示每条X射线投影在探测器上的位置坐标，R是圆轨道半径；图(b)是利用P-FDK方法把锥束投影重排成平行束的扇束投影，黑色实心点为X光源，在中心虚拟探测器上，重排后的投影数据由于同一角度下平行扇束到虚拟探测器的距离不相同，因此原平板探测器上一行投影对应到该虚拟探测器不再是水平的，而是在一段沿虚拟探测器中心行外凸的曲线上。以平板探测器为例，P-FDK的数据重排公式如公式(1)。本发明后续推导过程也皆以平板探测器为例，其他探测器类型例如柱面探测器等，都可以根据平板探测器做相应变形得到，不做赘述。

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,b)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},b)---\left(1\right)$

其中，θ表示经过P-FDK重排完后，其重排投影数据在角度方向的采样位置。(t，b)是经过重排后在中心虚拟探测器上的直角坐标系，表示重排后每条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-\mathrm{FDK}}(\mathrm{\θ},t,s)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{{\mathrm{sR}}^2}{R^2-t^2})---\left(2\right)$

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

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

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

发明内容

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

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

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

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

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

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

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

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{c\·R^2}{2(R^2-t^2)-R\sqrt{R^2-t^2}})$

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

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

${\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{{5R}^2-{4t}^2-4R\sqrt{R^2-t^2}+c^2}}\·P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)$

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

所述步骤对经过加权处理的投影数据沿水平方向进行滤波包括按下列公式进行滤波：

$g^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)={\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\⊗h\left(t\right)$

$=\underset{-c_0}{\overset{c_0}{\∫}}{\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\·h(t-t^{\′})d{t}^{\′}$

其中，g^C-FDK(θ，t，c)表示经过滤波的投影数据；表示卷积；h(t)表示滤波函数。

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

$f^{C-\mathrm{FDK}}(x,y,z)=\underset{0}{\overset{2\mathrm{\π}}{\∫}}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t(x,y,\mathrm{\θ}),c(x,y,z,\mathrm{\θ}))\mathrm{d\θ}$

其中，f^C-FDK(x，y，z)表示在x轴，y轴，z轴方向上的重建图像。相应地，本发明CT图像重建装置的技术方案包括：

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

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

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

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

所述重排单元按下列公式选取投影数据：

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{c\·R^2}{2(R^2-t^2)-R\sqrt{R^2-t^2}})$

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

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

${\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{{5R}^2-{4t}^2-4R\sqrt{R^2-t^2}+c^2}}\·P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)$

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

所述滤波单元按下列公式进行一维斜坡滤波：

$g^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)={\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\⊗h\left(t\right)$

$=\underset{-c_0}{\overset{c_0}{\∫}}{\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\·h(t-t^{\′})d{t}^{\′}$

其中，g^C-FDK(θ，t，c)表示经过滤波的投影数据；表示卷积；h(t)表示滤波函数。

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

$f^{C-\mathrm{FDK}}(x,y,z)=\underset{0}{\overset{2\mathrm{\π}}{\∫}}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t(x,y,\mathrm{\θ}),c(x,y,z,\mathrm{\θ}))\mathrm{d\θ}$

其中，f^C-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是图1(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(a)所示，首先，定义平面内锥束X射线源的扫描路径为一个圆周R是圆轨道半径；β表示光源点对应的角度参量；O是坐标原点和圆轨道的中心，也就是旋转中心。P(β，a，b)表示X射线照射被扫描物体后在平板探测器1上采集的投影数据，其中(a，b)表示二维平板探测器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的直线上。当然，OC曲线的曲率可以不准确等于扫描圆轨道的曲率，可以在圆轨道曲率的50％～200％范围之内，也即在此范围内的投影数据同样可以达到近似本发明的效果。

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

$c=\frac{\mathrm{SC}}{\mathrm{SP}}\·b=\frac{2\sqrt{R^2-t^2}-R}{R^2/\sqrt{R^2-t^2}}\·b=\frac{2(R^2-t^2)-R\sqrt{R^2-t^2}}{R^2}\·b---\left(3\right)$

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

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

${\mathrm{SO}}^{\′}=\sqrt{R^2-t^2}---\left(4\right)$

$\mathrm{SP}=R^2/\sqrt{R^2-t^2}---\left(5\right)$

$\mathrm{SC}=2\sqrt{R^2-t^2}-R---\left(6\right)$

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

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{c\·R^2}{2(R^2-t^2)-R\sqrt{R^2-t^2}})---\left(7\right)$

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

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

$\cos \mathrm{\γ}=\frac{\mathrm{SC}}{\sqrt{{\mathrm{SC}}^2+c^2}}=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{{5R}^2-{4t}^2-4R\sqrt{R^2-t^2}+c^2}}---\left(8\right)$

即对投影进行加权处理的公式为：

${\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{{5R}^2-{4t}^2-4R\sqrt{R^2-t^2}+c^2}}\·P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)---\left(9\right)$

然后，对完成上述处理的投影数据沿水平方向进行一维斜坡滤波：

$g^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)={\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\⊗h\left(t\right)$

$=\underset{-c_0}{\overset{c_0}{\∫}}{\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\·h(t-t^{\′})d{t}^{\′}---\left(10\right)$

上式中，表示卷积，而h(t)是滤波函数，一般采用Ramp滤波器。

最后，将滤波后的投影数据沿射线方向进行三维反投影即可重建出被扫描物体的三维CT图像，反投影公式如下：

$f^{C-\mathrm{FDK}}(x,y,z)=\underset{0}{\overset{2\mathrm{\π}}{\∫}}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t(x,y,\mathrm{\θ}),c(x,y,z,\mathrm{\θ}))\mathrm{d\θ}---\left(11\right)$

其中，探测器上的投影位置计算公式如下：

t(x，y，θ)＝ycosθ-xsinθ            (12)

$c(x,y,z,\mathrm{\θ})=\frac{z\·(2\sqrt{R^2-t^2}-R)}{\sqrt{R^2-t^2}+x\cos \mathrm{\θ}+y\sin \mathrm{\θ}}---\left(13\right)$

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

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

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

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

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

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

优选地，所述重排单元1按下列公式选取投影数据：

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{c\·R^2}{2(R^2-t^2)-R\sqrt{R^2-t^2}})$

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

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

${\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{{5R}^2-{4t}^2-4R\sqrt{R^2-t^2}+c^2}}\·P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)$

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

优选地，所述滤波单元3按如下公式进行滤波：

$g^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)={\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\⊗h\left(t\right)$

$=\underset{-c_0}{\overset{c_0}{\∫}}{\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\·h(t-t^{\′})d{t}^{\′}$

其中，g^C-FDK(θ，t，c)表示经过滤波的投影数据；表示卷积；h(t)表示滤波函数。

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

$f^{C-\mathrm{FDK}}(x,y,z)=\underset{0}{\overset{2\mathrm{\π}}{\∫}}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t(x,y,\mathrm{\θ}),c(x,y,z,\mathrm{\θ}))\mathrm{d\θ}$

其中，f^C-FDK(x，y，z)表示在x轴，y轴，z轴方向上的重建图像。

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

本申请的发明人利用三维Shepp-Logan头模型进行了数值模拟实验，并与FDK方法和T-FDK方法进行了实验对比来验证本发明的技术方案。在数值模拟实验中，三维头模型被限制在1mm的球体内，模型中心即为CT扫描的旋转中心，X光源到旋转中心距离为4mm，X光源到探测器的距离为8mm，平板探测器尺寸为4mm×4mm，探测单元数量为256×256，在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-FDK(θ，t，c)。而OC曲线则是和圆扫描轨道弧度相同的一段圆弧，此处OC曲线的选取不是唯一的，可以在该曲率附近做小范围的改变，也同样应该在本发明专利的权利要求范围内。

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

Claims (8)

1.一种CT图像重建方法，其特征在于，包括步骤：

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

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

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

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

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

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{c\·R^2}{2(R^2-t^2)-R\sqrt{R^2-t^2}})$

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

2.如权利要求1所述的CT图像重建方法，其特征在于，所述步骤对所选取的投影数据进行加权处理包括按下列公式对所选取的投影数据进行处理：

${\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{5R^2-4t^2-4R\sqrt{R^2-t^2}+c^2}}\·P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)$

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

3.如权利要求2所述的CT图像重建方法，其特征在于，所述步骤对经过加权处理的投影数据沿水平方向进行滤波包括按下列公式进行滤波：

$\begin{array}{c}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)={\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\⊗h\left(t\right)\\ =\underset{-c_0}{\overset{c_0}{\∫}}{\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\·h(t-t^{\′}){\mathrm{dt}}^{\′}\end{array}$

其中，g^C-FDK(θ,t,c)表示经过滤波的投影数据；表示卷积；h(t)表示滤波函数。

4.如权利要求3所述的CT图像重建方法，其特征在于，所述步骤对滤波后的投影数据沿射线方向进行三维反投影包括按下列公式进行三维反投影：

$f^{C-\mathrm{FDK}}(x,y,z)=\underset{0}{\overset{2\mathrm{\π}}{\∫}}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t(x,y,\mathrm{\θ}),c(x,y,z,\mathrm{\θ}))\mathrm{d\θ}$

其中，f^C-FDK(x,y,z)表示在x轴,y轴,z轴方向上的重建图像。

5.一种CT图像重建装置，其特征在于，包括：

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

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

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

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

所述重排单元按下列公式选取投影数据：

$P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=P(\mathrm{\θ}-\arcsin \frac{t}{R},\frac{\mathrm{tR}}{\sqrt{R^2-t^2}},\frac{c\·R^2}{2(R^2-t^2)-R\sqrt{R^2-t^2}})$

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

6.如权利要求5所述的CT图像重建装置，其特征在于，所述加权单元按下列公式对所选取的投影数据进行处理：

${\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)=\frac{2\sqrt{R^2-t^2}-R}{\sqrt{5R^2-4t^2-4R\sqrt{R^2-t^2}+c^2}}\·P^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)$

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

7.如权利要求6所述的CT图像重建装置，其特征在于，所述滤波单元按下列公式进行一维斜坡滤波：

$\begin{array}{c}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)={\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\⊗h\left(t\right)\\ =\underset{-c_0}{\overset{c_0}{\∫}}{\tilde{P}}^{C-\mathrm{FDK}}(\mathrm{\θ},t,c)\·h(t-t^{\′}){\mathrm{dt}}^{\′}\end{array}$

其中，g^C-FDK(θ,t,c)表示经过滤波的投影数据；表示卷积；h(t)表示滤波函数。

8.如权利要求7所述的CT图像重建装置，其特征在于，所述反投影单元按下列公式进行三维反投影：

$f^{C-\mathrm{FDK}}(x,y,z)=\underset{0}{\overset{2\mathrm{\π}}{\∫}}{g}^{C-\mathrm{FDK}}(\mathrm{\θ},t(x,y,\mathrm{\θ}),c(x,y,z,\mathrm{\θ}))\mathrm{d\θ}$

其中，f^C-FDK(x,y,z)表示在x轴,y轴,z轴方向上的重建图像。