CN105069823A - 基于非对称横向双边截断投影数据的扇束ct重建方法 - Google Patents

Abstract

本发明公开了一种基于非对称横向双边截断投影数据的扇束CT重建方法，该方法根据扇束投影数据的冗余性和对称性原理补充缺失的投影数据，使投影中心两侧的数据双边对称，消除投影中心偏置带来的影响；然后采用平滑延展法使两侧数据平滑过渡为零值，消除投影截断带来的影响；最后采用滤波反投影算法重建出断层图像。本发明基于扇束CT投影数据的冗余性和对称性原理，不需要对原始投影数据进行重排，一方面消除了投影数据非对称截断造成的重建图像伪影，另一方面充分利用了投影信息，扩大了重建区域。

Description

基于非对称横向双边截断投影数据的扇束CT重建方法

技术领域

本发明涉及扇束CT投影数据的重建方法，更特别地说，是指一种适用于投影数据横向双边截断且投影中心两侧数据非对称情况下的扇束CT重建方法，可用于医学和工业领域的二维计算机断层成像。

背景技术

扇束CT(Fan-beamComputedTomography，扇束计算机断层成像术)是一种先进的非接触、非破坏性二维CT检测手段，在医疗、航空、航天、船舶等领域有着广泛的应用。扇束CT扫描原理如图1所示，X射线源1发出的扇束射线2全面覆盖被扫描断层3，被扫描断层3绕旋转中心点4旋转360度；线阵探测器5采集不同旋转角度下穿过被扫描断层3的射线强度信号，该射线强度信号经过线阵探测器5的A/D转换模块形成数字信号并传至PC机，该数字信号被称为投影数据；PC机上的重建软件会对投影数据进行滤波反投影，重建出被扫描断层3的层析图像。

在对不同尺寸的样品实施检测过程中，受探测器尺寸或辐照范围的限制，检测样品的投影会超出探测器的有效成像区域，造成投影数据的缺失；另外，由于样品旋转中心点的投影很难做到位于探测器的中心点，造成了投影中心两边数据区域的不对称，以上两种因素共同作用导致了实际检测中产生了非对称的双边截断投影数据。这种非对称的双边截断投影数据会给重建结果带来严重的伪影，影响样品断层结构信息的准确判读。扇束CT非对称横向双边截断的扫描原理如图2所示。X射线源1发出扇束射线2，被扫描断层3绕旋转中心点4转动一周，旋转中心点4在线阵探测器5上的投影位置记为s₀。线阵探测器5的左边界点记为A点、右边界点记为B点，线阵探测器5的中心点记为s_中，所述s_中也是AB连线的中心点。由于被扫描断层的尺寸较大(记为大尺寸被扫描断层31)，而线阵探测器5的长度有限，导致扇束射线2不能对大尺寸被扫描断层31实现全覆盖，并且旋转中心点4在线阵探测器5上的投影位置s₀很难与线阵探测器5的中心点s_中重合，这两个因素共同作用造成了线阵探测器5无法采集到大尺寸被扫描断层31的完整投影数据，实际采集到的是双边截断的投影数据，即图2中61所指的曲线段。而且，投影数据左右截断的长度不等(即AC长度与BD长度不等)，因此这里将线阵探测器5实际采集到的投影数据61称为非对称双边截断投影。大尺寸被扫描断层31不能完全被线阵探测器5采集时，线阵探测器5左延展的扇束射线点记为C点，右延展的扇束射线点记为D点，AC连线的长度是投影数据左截断的长度，BD连线的长度是投影数据右截断的长度。

若要实现对大尺寸被扫描断层31的全覆盖，则需加大线阵探测器5的长度，使扇束射线2扩展为在图2中虚线所示，这里记为虚拟扇束射线21。虚拟扇束射线21对应的探测器称为虚拟探测器51，虚拟探测器51的左边界点即为线阵探测器5左延展的扇束射线点C、右边界点即为线阵探测器5右延展的扇束射线点D。虚拟探测器51采集到的投影数据为大尺寸被扫描断层31的完整投影信息，即图2中6所指的虚线段，称为完整投影数据。线阵探测器5实际采集到的投影数据61即为对完整投影数据6进行双边截断，并且左截断数据长度(AC连线长)不等于右截断数据长度(BD连线长度)。

在扇束CT扫描过程中，被扫描断层3每转过一个固定角度(该角度一般为0.1°～1°)，线阵探测器5采集一行投影数据，该行数据所在的坐标轴为线阵探测器坐标系轴8，记为s；被扫描断层3累计转过的角度称为投影角，记为θ，将投影角θ映射到线阵探测器坐标系轴9，构建一平面坐标系，如图3(a)、图3(b)所示。将所有投影角θ下线阵探测器5采集的行数据合成为一个二维投影矩阵CTP，则有：

$CTP=\left[\begin{array}{cccccc}p({\mathrm{\θ}}_1,s_1)& p({\mathrm{\θ}}_1,s_2)& \mathrm{...}& p({\mathrm{\θ}}_1,s_i)& \mathrm{...}& p({\mathrm{\θ}}_1,s_M)\\ p({\mathrm{\θ}}_2,s_1)& p({\mathrm{\θ}}_2,s_2)& \mathrm{...}& p({\mathrm{\θ}}_2,s_i)& \mathrm{...}& p({\mathrm{\θ}}_2,s_M)\\ .& .& & .& & .\\ .& .& & .& & .\\ .& .& & .& & .\\ p({\mathrm{\θ}}_j,s_1)& p({\mathrm{\θ}}_j,s_2)& \mathrm{...}& p({\mathrm{\θ}}_j,s_j)& \mathrm{...}& p({\mathrm{\θ}}_j,s_M)\\ .& .& & .& & .\\ .& .& & .& & .\\ .& .& & .& & .\\ p({\mathrm{\θ}}_N,s_1)& p({\mathrm{\θ}}_N,s_2)& \mathrm{...}& p({\mathrm{\θ}}_N,s_i)& \mathrm{...}& p({\mathrm{\θ}}_N,s_M)\end{array}\right],$ 该二维投影矩阵CTP称为非对称截断投影图，如3(a)、图3(b)中12所示。投影图12内任意一点的投影数据记为p(θ_j,s_i)，其中θ_j为投影角θ的任意一个取值，j＝1,2,3......N,N代表被扫描断层3在360°范围内的总的投影角数目，也就是投影图12的高度，一般取值为720～1800；s_i为线阵探测器坐标系轴s上的任意一坐标，也代表线阵探测器5上任意一探测单元的坐标值，i＝1,2,3......M，M代表探测单元的总数目，也就是投影图12的宽度。

图1所示的扫描模式下的图像重建就是基于二维投影矩阵CTP而实现。但是由于所述的CTP为非对称截断投影图，如果利用传统的滤波反投影重建方法，精确重建的区域为图2中71所示的区域，称为传统重建区域。而且，由于出现了投影数据的双边截断，层析图像的边缘会出现亮度不均的伪影，影响重构信息的准确判读。

发明内容

为了解决现有扇束CT系统出现双边截断投影数据不能完全重构层析图像的问题，本发明的创新点就是提出了新的重建方法，使得精确重建区域扩大为图2中72所示的区域。首先基于投影数据的冗余性和对称性原理补充缺失的投影数据，使投影中心两侧的数据双边对称，消除投影中心偏置带来的影响；然后采用平滑延展法使两侧数据平滑过渡为零值，消除投影截断带来的影响；最后采用滤波反投影算法重建出断层图像。利用本发明提出的方法，一方面消除了投影数据非对称截断造成的重建图像伪影；另一方面充分利用了投影信息，使得精确重建区域扩展至如图2中72所示的区域，该区域称为扩展重建区域。因此，本发明方法基于投影数据的冗余性和对称性原理，对非对称截断投影图进行对称延展，用以实现投影中心两侧数据的双边对称，消除投影中心偏置带来的影响。对称延展区任意一点的投影值等于该点在投影中心另一侧的共轭镜像点的投影值；然后采用平滑延展法使得对称延展投影图两侧数据平滑过渡为零值，消除投影截断带来的影响；最后采用滤波反投影重建得到被扫描样品的断层图像。

本发明基于非对称横向双边截断投影数据的扇束CT重建方法的优点在于：

①解决了被扫描样品大小受探测器尺寸限制的问题，实现了利用小尺寸探测器对大尺寸样品层析扫描的目标。

②基于扇束CT投影数据冗余性和对称性原理，不需要对原始投影数据进行重排，直接利用对称延展消除投影旋转中心的偏置，从而避免了重排中插值运算带来的数据精度损失。

③消除了由于投影数据截断造成的重建伪影，并且充分利用了投影数据的信息量，扩大了重建范围。

④本发明适用于任何具有旋转扫描功能的二维CT系统，且实现中无需对CT扫描系统进行硬件改造，可嵌入现有二维CT系统中作为辅助升级模块，有效提高了CT系统的检测能力。

附图说明

图1是扇束CT扫描原理图。

图2是非对称横向双边截断CT扫描原理图。

图3(a)是右对称延展示意图。

图3(b)是左对称延展示意图。

图4是平滑延展示意图。

图5(a)是对称延展后的投影图。

图5(b)是平滑延展后的投影图。

图6(a)是基于对称平滑延展投影的重建结果。

图6(b)是基于截断投影的重建结果。

1.X射线源	2.扇束射线	21.虚拟扇束射线
			3.被扫描断层	31.大尺寸被扫描断层	4.旋转中心点
5.线阵探测器	51.虚拟探测器	6.完整投影数据
			61.截断投影数据	71.传统重建区域	72.扩展重建区域
8.探测器坐标轴	9.投影角坐标轴	10.投影中心线
			11.右对称延展区	12.非对称截断投影图	13.左对称延展区
14.对称延展投影图	15.左平滑延展区	16.右平滑延展区

具体实施方式

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

如图3(a)所示，非对称截断投影图12的横坐标为探测器坐标系轴s，探测器由若干个独立的探测单元构成，某一个探测单元坐标记为s_i，i＝1,2,3......M，M代表探测单元的总数目，也就是非对称截断投影图的宽度，i为探测单元的标识号；纵坐标为投影角坐标轴θ，被扫描断层的某一个投影角记为θ_j，j＝1,2,3......N,N代表被扫描断层3在360°范围内的总的投影角数目，也就是非对称截断投影图的高度，一般取值为720～1800，j为投影角的标识号。由图1可知，所有投影角下旋转中心点4的投影位置固定为s₀，因此在非对称截断投影图12上旋转中心点投影的位置为过s₀平行于投影角坐标轴9的直线，该直线记为投影中心线10。结合图2，由于投影的非对称截断，导致投影图在探测器坐标系轴8上的中点s_中不与s₀重合。为了得到图2中72所示的拓展重建区域，就需将投影图12沿s轴方向进行延展，使投影图在探测器坐标系轴8上的中点位于s₀处，该操作定义为对称延展。当As₀＞s₀B时，对投影图进行右对称延展，右对称延展后的投影图记为CTP^R；当As₀＜s₀B时，对投影图进行左对称延展，左对称延展后的投影图记为CTP^L。当As₀＝s₀B时，则不需要对投影图进行延展。As₀表示A点与s₀点的连线的长度，s₀B表示s₀点与B点的连线的长度。

图3(a)为右对称延展示意图，由As₀＞s₀B可知，需要对投影图12的右边补充一定的数据，补充数据的区域称为右对称延展区11，该区域任意一点的投影值记为p^R(θ_j,s_t)，其中s_t∈[s_M+12s₀]，s_t为位于右对称延展区11内的任意一探测单元坐标。根据扇束CT扫描几何原理推导可知，在投影中心线10的左侧区域内可以找到与p^R(θ_j,s_t)等值的投影点，该点称为左侧区域共轭镜像点，记作p^R(θ_R-conj,s_R-conj)，θ_R-conj表示左侧区域内的投影角度的投影点，s_R-conj表示左侧区域内的行投影数据的投影点。因此只要找到该左侧区域共轭镜像点，将其值赋给p^R(θ_j,s_t)，即可完成对右对称延展区的数据填充。左侧区域共轭镜像点p^R(θ_R-conj,s_R-conj)的定位可通过θ_j和s_t求得，计算式如下：

${\mathrm{\θ}}_{R-conj}={\mathrm{\θ}}_j+\mathrm{\π}+2arctan\left(\frac{|s_t-s_0|}{D}\right)---\left(1\right)$

s_R-conj＝2s₀-s_t

D为射线源到探测器的距离，一般取值为800mm～1500mm，该值在CT系统出厂时给出。

右对称延展区11内的数据补充后，即可得到右对称延展后的投影图CTP^R，投影图CTP^R的宽度为2s₀，高度为N，投影中心线位于s₀处。右对称延展投影图CTP^R矩阵表达式如下：

${\mathrm{CTP}}^R=\left[\begin{array}{ccccccccccc}p({\mathrm{\θ}}_1,s_1)& p({\mathrm{\θ}}_1,s_2)& \mathrm{...}& p({\mathrm{\θ}}_1,s_i)& \mathrm{...}& p({\mathrm{\θ}}_1,s_M)& p^R({\mathrm{\θ}}_1,s_{M+1})& \mathrm{...}& p^R({\mathrm{\θ}}_1,s_t)& \mathrm{...}& p^R\left({\mathrm{\θ}}_1{s}_{2s_0}\right)\\ p({\mathrm{\θ}}_2,s_1)& p({\mathrm{\θ}}_2,s_2)& \mathrm{...}& p({\mathrm{\θ}}_2,s_i)& \mathrm{...}& p({\mathrm{\θ}}_2,s_M)& p^R({\mathrm{\θ}}_2,s_{M+1})& \mathrm{...}& p^R({\mathrm{\θ}}_2,s_t)& \mathrm{...}& p^R({\mathrm{\θ}}_2,s_{2s_0})\\ .& .& & .& & .& .& & .& & .\\ .& .& & .& & .& .& & .& & .\\ .& .& & .& & .& .& & .& & .\\ p({\mathrm{\θ}}_j,s_1)& p({\mathrm{\θ}}_j,s_2)& \mathrm{...}& p({\mathrm{\θ}}_j,s_i)& \mathrm{...}& p({\mathrm{\θ}}_j,s_M)& ({\mathrm{\θ}}_j,s_{M+1})& \mathrm{...}& p^R({\mathrm{\θ}}_j,s_t)& \mathrm{...}& p^R({\mathrm{\θ}}_j,s_{2s_0})\\ .& .& & .& & .& .& & .& & .\\ .& .& & .& & .& .& & .& & .\\ .& .& & .& & .& .& & .& & .\\ p({\mathrm{\θ}}_N,s_1)& p({\mathrm{\θ}}_N,s_2)& \mathrm{...}& p({\mathrm{\θ}}_N,s_i)& \mathrm{...}& p({\mathrm{\θ}}_N,s_M)& p^R({\mathrm{\θ}}_N,s_{M+1})& \mathrm{...}& p^R({\mathrm{\θ}}_N,s_t)& \mathrm{...}& p^R({\mathrm{\θ}}_N,s_{2s_0})\end{array}\right]$

同理，当As₀＜s₀B时，需要对投影图CTP进行左对称延展。图3(b)为左对称延展示意图。对投影图12的左边补充一定的数据，补充数据的区域称为左对称延展区13，该区域任意一点的投影值记为p^L(θ_j,s_k)，其中s_k∈[12s₀-s_M]，s_k为位于左对称延展区13内的任意一探测单元坐标。根据扇束CT扫描几何原理推导可知，在投影中心线10的右侧区域内可以找到与p^L(θ_j,s_k)等值的投影点，该点称为右侧区域共轭镜像点，记作p^L(θ_L-conj,s_L-conj)，θ_L-conj表示右侧区域内的投影角度的投影点，s_L-conj表示右侧区域内的行投影数据的投影点。因此只要找到该右侧区域共轭镜像点，将其值赋给p^L(θ_j,s_k)，即可完成对左对称延展区的数据填充。右侧区域共轭镜像点p^L(θ_L-conj,s_L-conj)的定位可通过θ_j和s_k求得，计算式如下：

${\mathrm{\θ}}_{L-conj}={\mathrm{\θ}}_j+\mathrm{\π}+2arctan\left(\frac{|s_k-s_0|}{D}\right)---\left(2\right)$

s_L-conj＝2s₀-s_k

左对称延展区13内的数据补充后，即可得到左对称延展后的投影图CTP^L。投影图CTP^L的宽度为2s₀，高度为N，投影中心线位于s₀处。左对称延展投影图CTP^L矩阵表达式如下：

${\mathrm{CTP}}^L=\left[\begin{array}{ccccccccccc}{p}^L({\mathrm{\θ}}_1,s_1)& \mathrm{...}& p^L({\mathrm{\θ}}_1,s_k)& \mathrm{...}& p^L({\mathrm{\θ}}_1,2s_0-s_M)& p({\mathrm{\θ}}_1,s_1)& p({\mathrm{\θ}}_1,s_2)& \mathrm{...}& p({\mathrm{\θ}}_1,s_i)& \mathrm{...}& p({\mathrm{\θ}}_1,s_M)\\ p^L({\mathrm{\θ}}_2,s_1)& \mathrm{...}& p^L({\mathrm{\θ}}_2,s_k)& \mathrm{...}& p^L({\mathrm{\θ}}_2,2s_0-s_M)& p({\mathrm{\θ}}_2,s_1)& p({\mathrm{\θ}}_2,s_2)& \mathrm{...}& p({\mathrm{\θ}}_2,s_i)& \mathrm{...}& p({\mathrm{\θ}}_2,s_M)\\ .& .& & .& .& .& & .& & .& \\ .& .& & .& .& .& & .& & .& \\ .& .& & .& .& .& & .& & .& \\ p^L({\mathrm{\θ}}_j,p_1)& \mathrm{...}& p^L({\mathrm{\θ}}_j,s_k)& \mathrm{...}& p^L({\mathrm{\θ}}_j,2s_0-s_M)& p({\mathrm{\θ}}_j,s_1)& p({\mathrm{\θ}}_j,s_2)& \mathrm{...}& p({\mathrm{\θ}}_j,s_i)& \mathrm{...}& p({\mathrm{\θ}}_j,s_M)\\ .& .& & .& .& .& & .& & .& \\ .& .& & .& .& .& & .& & .& \\ .& .& & .& .& .& & .& & .& \\ p^L({\mathrm{\θ}}_N,s_1)& \mathrm{...}& p^L({\mathrm{\θ}}_N,s_k)& \mathrm{...}& p^L({\mathrm{\θ}}_N,2s_0-s_M)& p({\mathrm{\θ}}_N,s_1)& p({\mathrm{\θ}}_N,s_2)& \mathrm{...}& p({\mathrm{\θ}}_N,s_i)& \mathrm{...}& p({\mathrm{\θ}}_N,s_M)\end{array}\right]$

然而，对称延展后的投影图(CTP^R或CTP^L)依旧为截断投影数据，即该投影图的第一列和最后一列的投影值不为零，如果直接利用该截断投影数据进行滤波反投影重建，则会在重建图像中产生伪影。本发明提出了改进方法，在对称延展的基础上再进行一次延展，对CTP^R或CTP^L两边填充数据，使得两边的投影数据平滑过渡到零，该步操作称为平滑延展。如图4所示，对称延展后的投影图(CTP^R或CTP^L)记为14，填充数据的区域分为左右两部分，分别称为左平滑延展区15和右平滑延展区16。左右平滑延展区的宽度均为L，高度均为N。L一般取值范围为[M/8M/4]。

左平滑延展区15内的任意一个位置点的投影值记为p^left(θ_j,s_b)，b＝1,2,......L。则有对于右对称延展投影图，左平滑延展区15内的投影值为：

$p^{left}({\mathrm{\θ}}_j,s_b)=\max \{\sin \frac{(s_b-1)\mathrm{\π}}{2(s_L-1)}\×\[2p({\mathrm{\θ}}_j,s_1)-p({\mathrm{\θ}}_j,s_L-s_b+2)\],0\}---\left(3\right)$

在本发明中，对于左对称延展投影图，左平滑延展区15内的投影值为：

$p^{left}({\mathrm{\θ}}_j,s_b)=max\{sin\frac{(s_b-1)\mathrm{\π}}{2(s_L-1)}\×\[2p^L({\mathrm{\θ}}_j,s_1)-p^L({\mathrm{\θ}}_j,s_L-s_b+2)\],0\}---\left(4\right)$

s_b为s在探测器坐标系轴8上[0s_L]区域内的取值。

右平滑延展区16内的任意一个位置的投影值记为p^right(θ_j,s_c)，c＝1,2.......L。则有对于右对称延展投影图，右平滑延展区16内的投影值为：

$p^{right}({\mathrm{\θ}}_j,s_c)=max\{cos\frac{(s_c-1)\mathrm{\π}}{2(s_L-1)}\×\[2p^R({\mathrm{\θ}}_j,2s_0)-p^R({\mathrm{\θ}}_j,2s_0-s_c)\],0\}---\left(5\right)$

在本发明中，对于左对称延展投影图，右平滑延展区16内的投影值为：

$p^{right}({\mathrm{\θ}}_j,s_c)=\max \{\cos \frac{(s_c-1)\mathrm{\π}}{2(s_L-1)}\×\[2p({\mathrm{\θ}}_j,s_M)-p({\mathrm{\θ}}_j,s_M-s_c)\],0\}---\left(6\right)$

s_c为s在探测器坐标系轴8上[2s₀+s_L2s₀+2s_L]区域内的取值。

在本发明中，通过左平滑延展区15及右平滑延展区16各自得到的投影图构建得到平滑延展后的投影图，记为CTP^smooth。利用所述CTP^smooth进行滤波反投影重建，即可重构出图2中的扩展重建区域72里的信息。

本发明是一种基于非对称横向双边截断投影数据的扇束CT重建方法，结合上述推导，该测量方法包括有下列实施步骤：

步骤一：采集扫描断层原始投影数据；

启动CT扫描系统，使得被扫描断层旋转360°，探测器获取被扫描断层的原始投影图，记为CTP^raw，所述CTP^raw是宽度为M、高度为N的二维投影矩阵；并记录旋转中心在探测器上的投影坐标s₀。

在本发明中，原始投影图CTP^raw的二维投影矩阵形式为：

p^raw(θ_j,s_i)表示原始投影图CTP^raw的任意一个点的投影值。

步骤二：获取参考投影图；

在CT扫描系统的扫描平台上移去被扫描物体，即射线源与探测器之间不放置任何物体，探测器采集一幅参考投影图，记为CTP^ref；显然，CTP^ref亦为宽度为M、高度为N的二维投影矩阵。利用参考投影图CTP^ref计算出射线的初始强度值I₀，计算方法为：

$I_0=\frac{1}{M\×N}\underset{j=1}{\overset{N}{\Σ}}\underset{i=1}{\overset{M}{\Σ}}{p}^0({\mathrm{\θ}}_j,s_i)---\left(7\right)$

p⁰(θ_j,s_i)表示考投影图CTP^ref的任意一个点的投影值。

在本发明中，参考投影图CTP^ref的二维投影矩阵形式为：

$\left[\begin{array}{cccccc}{p}^0({\mathrm{\θ}}_1,s_1)& p^0({\mathrm{\θ}}_1,s_2)& \mathrm{...}& p^0({\mathrm{\θ}}_1,s_i)& \mathrm{...}& p^0({\mathrm{\θ}}_1,s_M)\\ p^0({\mathrm{\θ}}_2,s_1)& p^0({\mathrm{\θ}}_2,s_2)& \mathrm{...}& p^0({\mathrm{\θ}}_2,s_i)& \mathrm{...}& p^0({\mathrm{\θ}}_2,s_M)\\ .& .& & .& & .\\ .& .& & .& & .\\ .& .& & .& & .\\ p^0({\mathrm{\θ}}_j,s_1)& p^0({\mathrm{\θ}}_j,s_2)& \mathrm{...}& p^0({\mathrm{\θ}}_j,s_i)& \mathrm{...}& p^0({\mathrm{\θ}}_j,s_M)\\ .& .& & .& & .\\ .& .& & .& & .\\ .& .& & .& & .\\ p^0({\mathrm{\θ}}_N,s_1)& p^0({\mathrm{\θ}}_N,s_2)& \mathrm{...}& p^0({\mathrm{\θ}}_N,s_i)& \mathrm{...}& p^0({\mathrm{\θ}}_N,s_M)\end{array}\right].$

步骤三：采用对数变换方法获取非对称截断投影图；

对原始投影图CTP^raw进行对数变换得到非对称截断投影图CTP，CTP的任意一点投影值p(θ_j,s_i)由下式计算：

$p({\mathrm{\θ}}_j,s_i)=\ln \frac{I_0}{p^{raw}({\mathrm{\θ}}_j,s_i)},j=1,2,\mathrm{3......}N,i=1,2,\mathrm{3......}M---\left(8\right)$

步骤四：对称延展处理；

根据s₀与非对称截断投影图CTP的宽度(探测器探测单元的总数目)M的大小关系来判断是否对称延展及延展区的位置，即：如果则进行右对称延展；

如果则进行左对称延展；如果则不进行对称延展。

步骤41：右对称延展处理；

当需要右对称延展时，补充数据的区域称为右对称延展区11，该区域任意一点的投影值记为p^R(θ_j,s_t)，其中s_t∈[s_M+12s₀]。根据扇束CT扫描几何原理推导可知，在投影中心线10的左侧区域内可以找到与p^R(θ_j,s_t)等值的投影点，该点称为左侧区域共轭镜像点，记作p^R(θ_R-conj,s_R-conj)。因此只要找到该左侧区域共轭镜像点，将其值赋给p^R(θ_j,s_t)即可完成对右对称延展区的数据填充。左侧区域共轭镜像点p^R(θ_R-conj,s_R-conj)的定位可通过θ_j和s_t求得，计算式如下：

${\mathrm{\θ}}_{R-conj}={\mathrm{\θ}}_j+\mathrm{\π}+2arctan\left(\frac{|s_t-s_0|}{D}\right)---\left(1\right)$

s_R-conj＝2s₀-s_t

步骤42：左对称延展处理；

当需要左对称延展时，补充数据的区域称为左对称延展区13，该区域任意一点的投影值记为p^L(θ_j,s_k)，其中s_k∈[12s₀-s_M]。根据扇束CT扫描几何原理推导可知，在投影中心线10的右侧区域内可以找到与p^L(θ_j,s_k)等值的投影点，该点称为右侧区域共轭镜像点，记作p^L(θ_L-conj,s_L-conj)。因此只要找到该右侧区域共轭镜像点，将其值赋给p^L(θ_j,s_k)即可完成对左对称延展区的数据填充。右侧区域共轭镜像点p^L(θ_L-conj,s_L-conj)的定位可通过θ_j和s_k求得，计算式如下：

${\mathrm{\θ}}_{L-conj}={\mathrm{\θ}}_j+\mathrm{\π}+2arctan\left(\frac{|s_k-s_0|}{D}\right)---\left(2\right)$

s_L-conj＝2s₀-s_k

D为射线源到探测器的距离，一般取值为800mm～1500mm，该值在CT系统出厂时给出。

步骤五：平滑延展处理；

对投影图CTP^R或CTP^L进行平滑延展。即在CTP^R或CTP^L的左右两边分别补充宽度为L，高度为N的扩展区，分别称为左平滑延展区15和右平滑延展区16。L一般取值范围为[M/8M/4]。

左平滑延展区15内的任意一个位置的投影值记为p^left(θ_j,s_b)，b＝1,2,......L。则有对于右对称延展投影图，左平滑延展区15内的投影值为：

$p^{left}({\mathrm{\θ}}_j,s_b)=max\{sin\frac{(s_b-1)\mathrm{\π}}{2(s_L-1)}\×\[2p({\mathrm{\θ}}_j,s_1)-p({\mathrm{\θ}}_j,s_L-s_b+2)\],0\}---\left(3\right)$

在本发明中，对于左对称延展投影图，左平滑延展区15内的投影值为：

$p^{left}({\mathrm{\θ}}_j,s_b)=max\{sin\frac{(s_b-1)\mathrm{\π}}{2(s_L-1)}\×\[2p^L({\mathrm{\θ}}_j,s_1)-p^L({\mathrm{\θ}}_j,s_L-s_b+2)\],0\}---\left(4\right)$

s_b为s在探测器坐标系轴8上[0s_L]区域内的取值。

右平滑延展区16内的任意一个位置的投影值记为p^right(θ_j,s_c)，c＝1,2.......L。则有对于右对称延展投影图，右平滑延展区16内的投影值为：

$p^{right}({\mathrm{\θ}}_j,s_c)=\max \{\cos \frac{(s_c-1)\mathrm{\π}}{2(s_L-1)}\×\[2p^R({\mathrm{\θ}}_j,2s_0)-p^R({\mathrm{\θ}}_j,2s_0-s_c)\],0\}---\left(5\right)$

在本发明中，对于左对称延展投影图，右平滑延展区16内的投影值为：

$p^{right}({\mathrm{\θ}}_j,s_c)=max\{cos\frac{(s_c-1)\mathrm{\π}}{2(s_L-1)}\×\[2p({\mathrm{\θ}}_j,s_M)-p({\mathrm{\θ}}_j,s_M-s_c)\],0\}---\left(6\right)$

s_c为s在探测器坐标系轴8上[2s₀+s_L2s₀+2s_L]区域内的取值。

在本发明中，通过左平滑延展区15及右平滑延展区16各自得到的投影图构建得到平滑延展后的投影图，记为CTP^smooth。利用所述CTP^smooth进行滤波反投影重建，即可重构出图2中的扩展重建区域72里的信息。

实施案例

(一)实验采用的扫描装置参数如下：

(1)射线源：日本产L9181S型130kV微焦点射线源，焦点尺寸5μm；

(2)线阵探测器：探测器像元数目为359，像元尺寸为0.2mm；

(3)射线源焦点与探测器的距离为696mm。

(二)实验步骤：

(A)将一被检测样品放置在CT系统的扫描平台上，启动扫描平台旋转带动被检测样品步进旋转360°，步进角为0.5°，共采集720行投影数据。将所有投影角度下探测器采集的行数据合成一个二维投影矩阵，该合成图像即为原始投影图CTP^raw。

(B)移去检测样品，在同样条件下获得参考图像CTP^ref，依据公式(7)计算出射线的初始强度值I₀，该案例中的I₀＝3620。根据公式(8)对CTP^raw进行对数变换得到非对称截断投影图CTP，测得s₀＝216.5。

(C)根据步骤四在投影图右侧补充数据，对称延展区域的宽度为75，由公式(1)求得称延展区域内的填充数据。右对称延展后的投影图CTP^R如图5(a)所示。

(D)根据步骤五对称延展后的投影图CTP^R进行平滑延展，平滑延展区域的宽度L＝M/6＝359/6≈60，并按照步骤六加权处理。平滑延展后的投影图CTP^smooth如图4(b)所示。

(E)对CTP^smooth进行滤波反投影重建，结果如图6(a)所示。图6(b)为利用原始的非对称双边截断投影图CTP的重建结果。对比图6(a)和图6(b)可以看出，本发明的方法有效扩大了重建区域，而且重建图像中没有伪影，更为清晰。

Claims (4)

1.一种基于非对称横向双边截断投影数据的扇束CT重建方法，所述CT扫描装置上至少包括有射线源(1)、扇束射线(2)和线阵探测器(5)，其特征在于：基于非对称横向双边截断投影数据的扇束CT重建有下列步骤；

步骤一：采集扫描断层原始投影数据；

启动CT扫描系统，使得被扫描断层旋转360°，探测器获取被扫描断层的原始投影图，记为CTP^raw，所述CTP^raw是宽度为M、高度为N的二维投影矩阵；并记录旋转中心在探测器上的投影坐标s₀；

步骤二：获取参考投影图；

在CT扫描系统的平台上移去被扫描物体，即射线源与探测器之间不放置任何物体，探测器采集一幅参考投影图，记为CTP^ref；显然，CTP^ref亦为宽度为M、高度为N的二维投影矩阵；利用参考投影图CTP^ref计算出射线的初始强度值I₀，计算方法为：

$I_0=\frac{1}{M\×N}\underset{j=1}{\overset{N}{\Σ}}\underset{i=1}{\overset{M}{\Σ}}{p}^0({\mathrm{\θ}}_j,s_i)---\left(7\right)$

步骤三：采用对数变换方法获取非对称截断投影图；

对CTP^raw进行对数变换得到非对称截断投影图CTP，CTP的任意一点投影值p(θ_j,s_i)由下式计算：

$p({\mathrm{\θ}}_j,s_i)=\ln \frac{I_0}{p^{raw}({\mathrm{\θ}}_j,s_i)},j=1,2,\mathrm{3......}N,i=1,2,\mathrm{3......}M---\left(8\right)$

p^raw(θ_j,s_i)表示原始投影图CTP^raw的任意一个点的投影值；

步骤四：对称延展处理；

根据s₀与投影图CTP的宽度M的大小关系来判断是否对称延展及延展区的位置，即：如果则进行右对称延展；如果则进行左对称延展；如果则不进行对称延展；

步骤41：右对称延展处理；

当需要右对称延展时，补充数据的区域称为右对称延展区(11)，该区域任意一点的投影值记为p^R(θ_j,s_t)，其中s_t∈[s_M+12s₀]；根据扇束CT扫描几何原理推导可知，在投影中心线(10)的左侧区域内可以找到与p^R(θ_j,s_t)等值的投影点，该点称为左侧区域共轭镜像点，记作p^R(θ_R-conj,s_R-conj)；因此只要找到该左侧区域共轭镜像点，将其值赋给p^R(θ_j,s_t)即可完成对右对称延展区的数据填充；左侧区域共轭镜像点p^R(θ_R-conj,s_R-conj)的定位可通过θ_j和s_t求得，计算式如下：

${\mathrm{\θ}}_{R-conj}={\mathrm{\θ}}_j+\mathrm{\π}+2arctan\left(\frac{|s_t-s_0|}{D}\right)---\left(1\right)$

s_R-conj＝2s₀-s_t

步骤42：左对称延展处理；

当需要左对称延展时，补充数据的区域称为左对称延展区(13)，该区域任意一点的投影值记为p^L(θ_j,s_k)，其中s_k∈[12s₀-s_M]；根据扇束CT扫描几何原理推导可知，在投影中心线10的右侧区域内可以找到与p^L(θ_j,s_k)等值的投影点，该点称为右侧区域共轭镜像点，记作p^L(θ_L-conj,s_L-conj)；因此只要找到该右侧区域共轭镜像点，将其值赋给p^L(θ_j,s_k)即可完成对左对称延展区的数据填充；右侧区域共轭镜像点p^L(θ_L-conj,s_L-conj)的定位可通过θ_j和s_k求得，计算式如下：

${\mathrm{\θ}}_{L-conj}={\mathrm{\θ}}_j+\mathrm{\π}+2arctan\left(\frac{|s_k-s_0|}{D}\right)---\left(2\right)$

s_L-conj＝2s₀-s_k

D为射线源到探测器的距离，一般取值为800mm～1500mm，该值在CT系统出厂时给出；

步骤五：平滑延展处理；

对投影图CTP^R或CTP^L进行平滑延展；即在CTP^R或CTP^L的左右两边分别补充宽度为L，高度为N的扩展区，分别称为左平滑延展区(15)和右平滑延展区(16)；L一般取值范围为[M/8M/4]；

左平滑延展区(15)内的任意一个位置的投影值记为p^left(θ_j,s_b)，b＝1,2,......L；则对于右对称延展投影图，左平滑延展区(15)内的投影值为：

$p^{left}({\mathrm{\θ}}_j,s_b)=max\{sin\frac{(s_b-1)\mathrm{\π}}{2(s_L-1)}\×\[2p({\mathrm{\θ}}_j,s_1)-p({\mathrm{\θ}}_j,s_L-s_b+2)\],0\}---\left(3\right)$

对于左对称延展投影图，左平滑延展区(15)内的投影值为：

$p^{left}({\mathrm{\θ}}_j,s_b)=max\{sin\frac{(s_b-1)\mathrm{\π}}{2(s_L-1)}\×\[2p^L({\mathrm{\θ}}_j,s_1)-p^L({\mathrm{\θ}}_j,s_L-s_b+2)\],0\}---\left(4\right)$

s_b为s在探测器坐标系轴(8)上[0s_L]区域内的取值；

右平滑延展区(16)内的任意一个位置的投影值记为p^right(θ_j,s_c)，c＝1,2.......L；则对于右对称延展投影图，右平滑延展区(16)内的投影值为：

$p^{right}({\mathrm{\θ}}_j,s_c)=max\{cos\frac{(s_c-1)\mathrm{\π}}{2(s_L-1)}\×\[2p^R({\mathrm{\θ}}_j,2s_0)-p^R({\mathrm{\θ}}_j,2s_0-s_c)\],0\}---\left(5\right)$

对于左对称延展投影图，右平滑延展区(16)内的投影值为：

$p^{right}({\mathrm{\θ}}_j,s_c)=\max \{\cos \frac{(s_c-1)\mathrm{\π}}{2(s_L-1)}\×\[2p({\mathrm{\θ}}_j,s_M)-p({\mathrm{\θ}}_j,s_M-s_c)\],0\}---\left(6\right)$

s_c为s在探测器坐标系轴(8)上[2s₀+s_L2s₀+2s_L]区域内的取值；

通过左平滑延展区(15)及右平滑延展区(16)各自得到的投影图构建得到平滑延展后的投影图CTP^smooth，利用所述CTP^smooth进行滤波反投影重建，即可重构出扩展重建区域(72)里的信息。

2.根据权利要求1所述的基于非对称横向双边截断投影数据的扇束CT重建方法，其特征在于：被检测样品步进旋转360度，步进角一般取0.1°～1°；射线源到探测器的距离为800mm～1500mm。

3.根据权利要求1所述的基于非对称横向双边截断投影数据的扇束CT重建方法，其特征在于：扇束射线不能全部覆盖被扫描断层，造成投影数据的双边截断。

4.根据权利要求1所述的基于非对称横向双边截断投影数据的扇束CT重建方法，其特征在于：被扫描断层旋转中心在探测器上的投影点不位于线阵探测器的中心，投影中心两侧的数据非对称。