CN101404088A - Ct图像重建的方法及系统 - Google Patents

Abstract

本发明公开了一种CT图像重建的方法，包括识别图像的有效数据边界，获得所述有效数据边界点行(列)号的最小值与最大值，对所述图像进行数据加权，将扇形束校正为对应的平行束，并通过卷积滤波消除所述图像中的伪影，利用所述行(列)号的最小值与最大值进行反投影重建。本发明还公开了一种CT图像重建的装置，包括边界识别模块、数据加权模块、滤波模块和重建模块。本发明通过在CT图像重建之前进行识别有效数据边界，之后仅对边界内的有效数据进行重建，在保证重建图像质量不变的前提下，使重建速度得到了提高。

Description

CT图像重建的方法及系统

技术领域

本发明涉及数字图像处理及模式识别领域，特别是CT图像重建的方法及系统。

背景技术

CT(Computed Tomography，计算机断层摄影术)作为一种非接触式无损探测技术，近年来在医学、工业无损探测和相关军事领域得到了广泛的应用。通过CT断层或3D重建，我们可以得到目标物体的大小、形状、尺寸等多方面信息。而实时CT断层或3D图像重建在医学领域如颅内手术或心脏手术中，以及相关军事领域中将占据越来越重要的位置。除了扫描速度，另一个影响实时CT图像重建的关键因素是重建速度。

目前，在CT的2D和3D重建领域应用最广泛的重建算法为FBP算法与FDK算法。FBP与FDK算法大致可分为3步：数据加权、卷积滤波和反投影重建。虽然这两种算法经过了一定程度的优化，但是由于它们进行反投影计算时浪费了很大一部分时间在背景数据上，而一般情况下使用者并不关心背景的重建，导致这两种算法的重建速度相对较慢。

发明内容

有鉴于此，本发明的目的在于提供CT图像重建的方法及系统，用于提高CT图像的重建速度。

为实现上述目的，本发明提供了一种CT图像重建的方法，包括：

识别图像的有效数据边界，获得所述有效数据边界点行(列)号的最小值与最大值，对所述图像进行数据加权，将扇形束校正为对应的平行束，并通过卷积滤波消除所述图像中的伪影，利用所述行(列)号的最小值与最大值进行反投影重建。

本发明还提供了一种CT图像重建的装置，包括：

边界识别模块，用于识别图像的有效数据边界，获得所述有效数据边界点行(列)号的最小值与最大值；

数据加权模块，用于对所述图像进行数据加权，将扇形束校正为对应的平行束；

滤波模块，用于通过卷积滤波消除所述图像中的伪影；

重建模块，用于利用所述行(列)号的最小值与最大值进行反投影重建。

本发明通过在CT图像3步重建之前增加了识别有效数据边界的步骤，即作为一种对投影数据的预处理，之后仅对边界内的有效数据进行重建，在保证重建图像质量不变的前提下，使重建速度得到了提高。通过理论分析和实验验证，重建速度具体提高的倍数为探测器尺寸与物体在探测器上投影最宽处尺寸之比的二次方。

附图说明

图1为本发明的实施例1中CT图像重建的方法流程图；

图2为本发明的实施例中投影坐标系和实物坐标系关系的示意图；

图3为本发明的实施例1中扇形束CT扫描物体的原理图；

图4为本发明的实施例2中海螺的3D投影数据在xy平面上的一张切片图；

图5为本发明实施例中两种重建算法在海螺重建质量上的直观对比图；

图6为现有技术提供的与本发明实施例提供的重建方法的重建质量对比示意图；

图7为采用本发明实施例的图像重建方法重建锥体的示意图；

图8为本发明实施例的图像重建方法对重建速度提高的效果示意图。

具体实施方式

本发明的实施例通过在CT图像3步重建之前增加了识别有效数据边界的步骤，即作为一种对投影数据的预处理，之后仅对边界内的有效数据进行重建，在保证重建图像质量不变的前提下，使重建速度得到了提高。

通过对FBP与FDK算法的研究表明，反投影重建部分所消耗的时间占整个重建时间的98％以上。可见，提高CT重建速度最有效的方式就是减少反投影重建部分所消耗的时间，可以通过减少反投影数据量来实现这一目标。

为了满足大部分需求，CT的探测器尺寸一般要比目标物体的尺寸大得多。这样，探测器就采集到了大量的无效数据，如背景数据，非目标物体数据等。如果可以仅对目标物体数据进行反投影重建运算，那么就可以节约大量的重建时间。本发明实施例在对图像进行重建之前，运用数字图像处理或模式识别的手段，如对投影数据二值化或边缘检测，找到二值化或边缘检测图像中像素值突变的点，记录其行(列)号的最大值与最小值作为数据边界。重建时仅对上述边界内的有效数据进行反投影运算，从而可以达到加快重建速度的目的。

在CT的反投影重建部分中，以2D重建FBP算法为例。设投影角度数为L，探测器尺寸为M，X射线穿过目标物体在探测器上投影为半径为N的圆形物体。在原始算法中，需要进行L×M²次循环运算。但是，其中的背景数据对重建结果并没有贡献。而如果只对目标物体数据进行反投影重建运算，那么循环次数将减少为L×N²。这种情况下，反投影重建部分的运算时间将减少为原来的(N/M)²。

用Y表示只对目标物体数据进行反投影重建运算的重建时间与原始算法重建时间之比，用X表示目标物体最宽处尺寸与探测器尺寸之比，则Y、X满足以下关系：

Y＝X²

本发明实施例就采用了这种只对目标物体数据进行重建运算的方法，通过以上公式能够得出本发明实施例理论上的提速效果。

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明作进一步的详细描述。

如图1所示，为本发明的实施例1中CT图像重建的方法流程图，该方法包括以下步骤：

步骤101、有效数据边界识别。识别有效数据的边界，仅对边界内的有效数据进行重建，则可以大大提高重建速度。首先读取图像数据，运用数字图像处理的方法处理图像，如二值化或边缘检测处理，寻找有效数据的边界，即像素值突变的点，保存所有突变点的行(列)号的最小值与最大值，该像素值突变的点即为有效数据的边界。当然不仅限于使用这两种方式进行图像处理，但上述两种是最直接的，并且运行效率较高。保存并返回行(列)号的最小值与最大值，可以让计算机只处理行(列)号的最小值与最大值之间的数据(即有效数据)，从而提高重建速度。

运用二值化的方法处理图像，这里以Otsu二值化方法为例：

该方法中，二值化最优阈值k′由下面的公式定义： ${\mathrm{\σ}}_B^2\left(k^{\′}\right)=\underset{1\≤k\≤L}{\max }{\mathrm{\σ}}_B^2\left(k\right),$ 这种方法通过计算类方差σ_B ²来确定k′。类方差σ_B ²可由下列公式算出： ${\mathrm{\σ}}_B^2={\stackrel{\‾}{\mathrm{\ω}}}_0{({\mathrm{\μ}}_0-{\mathrm{\μ}}_T)}^2+{\stackrel{\‾}{\mathrm{\ω}}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2={\stackrel{\‾}{\mathrm{\ω}}}_0{\stackrel{\‾}{\mathrm{\ω}}}_1{({\mathrm{\μ}}_0-{\mathrm{\μ}}_1)}^2.$ 其中ω₀、ω₁是类0与类1中的像素点的个数占总像素点个数的百分比；μ₀、μ₁是两类中像素值的平均数；μ_T为整幅图像总的像素值的平均数。其中，类0和类1由k划分，像素值小于k的像点归入类0，大于k的归入类1。

在MATLAB中，有算Otsu二值化最优阈值的函数：graythresh()。在VC中，亦有相应的开发包计算Otsu二值化最优阈值。

然后，利用上面计算的Otsu二值化最优阈值来对图像进行二值化处理。遍历整幅图像中每个像点，如果当前像点的像素值大于或等于Otsu二值化最优阈值，则将此点像素值置1，反之则置0。将行(列)号相同像点的像素值累加，组成一个一维矩阵，则矩阵中第一个和最后一个值为非0的点为像素值突变点。此时，背景数据(无效数据，非目标物体数据)像素值变为0，目标物体数据(有效数据)像素值变为1，则像素值突变点的行(列)号的最小值与最大值即为有效数据边界。

步骤102、数据加权。

步骤103、卷积滤波。根据一定的滤波函数对投影数据进行卷积运算，来消除反投影重建中可能出现的伪影。一个好的滤波函数对重建图像的质量有很大改善，本本发明实施例中在比较两种算法重建质量时使用了同一种滤波函数。

步骤104、反投影重建。对经过上面3步处理的数据进行Radon反变换来重建目标物体。Radon反变换的作用就是将图2旋转坐标系(即投影坐标系)XOY中的数据反变换到静止坐标系(即实物坐标系)xOy中，从而重建物体图像。以二维Radon反变换为例，它是用Radon反变换算子R^-1处理投影数据，使得[R^-1Pe’](β，S(β，x，y))＝μ(x，y)。Radon反变换算子R^-1可表示为R^-1＝1/2π[BHD]，其中D两个极变函数关于第一变量的偏导数算子；H为两个极变量函数关于第一变量的希尔伯特算子；B为两个变量的反投影算子；1/2π为归一化算子。

下面以扇形束FBP重建算法为例，对数据加权、卷积滤波和反投影重建的方法进行阐述。

如图2所示，xOy为实物坐标系，XOY为投影坐标系。实物坐标系位置固定不变，投影坐标系随射线源和探测器的绕O点旋转而旋转，如图3所示，其中，射线源A在Y轴上，探测器所在X轴上。CT图像重建要做的工作就是把XOY坐标系上获得的投影图像通过Radon反变换，将其还原到实物坐标系xOy上，从而重建物体图像。

通过图2，可以得出两坐标系的变换关系为：

X＝xcosβ+ysinβ

Y＝-xsinβ+ycosβ

其中，β为两坐标系间夹角，也称投影角度。

在图3的等距扇束投影图中，A为射线源，其绝对位置由β确定；AP为某一射线，其绝对位置由(β，S)确定，S为X轴与该射线交点的距离，即AQ，当β已经确定的时候，S唯一确定一条射线；而射线AP上的一点P，由(β，S，U)三个变量确定。U为该点到Y轴垂线的交点到射线发射点的距离，即AC，当β和S已经确定时，U唯一确定这条射线上的一个点。

则重建公式如下：

$\mathrm{\μ}(x,y)={\∫}_0^{2\mathrm{\π}}\frac{D^2}{U{(\mathrm{\β},x,y)}^2}{P}^{e^{\′}}(\mathrm{\β},S(\mathrm{\β},x,y))\mathrm{d\β}$

其中，μ(x，y)为物体在实物坐标系中的密度分布函数，即代表真实的物体图像。

AO为射线源到探测器距离，令AO＝D，则有：

U(β，x，y)＝D+(-Y)＝D+xsinβ-ycosβ；

$S(\mathrm{\β},x,y)=D\tan \mathrm{\δ}=D\frac{X}{U}=\frac{D(x\cos \mathrm{\β}+y\sin \mathrm{\β})}{D+x\sin \mathrm{\β}-y\cos \mathrm{\β}};$

$P^e(\mathrm{\β},S)=P(\mathrm{\β},S)\cos \mathrm{\δ}=P(\mathrm{\β},S)\frac{D}{\sqrt{D^2+S^2}};$

P^e′(β，S)＝P^e(β，S)*g(S)

其中，g(S)为已知滤波函数。

在整个重建过程中，由P计算Pe的过程为数据加权，目的是在扇束或锥束扫描中将非中心射线校正为与中心射线平行的射线，从而校正由其产生的投影，使其可按平行束滤波反投影算法进行图像重建；由Pe计算Pe’的过程为卷积滤波，目的是为了消除CT重建图像中的伪影，主要影响重建图像的质量；而计算U(β，x，y)，S(β，x，y)，并由D，S(β，x，y)和Pe’(β，S(β，x，y))来积分计算出μ(x，y)的整个过程称为反投影重建。至此，就得到了完整的物体图像，完成CT图像重建。

在以上实施例中，通过识别有效数据的边界，之后仅对边界内的有效数据进行重建，在保证重建图像质量不变的前提下，使重建速度得到了提高。通过实验验证，重建速度具体提高的倍数为探测器尺寸与物体在探测器上投影最宽处尺寸之比的二次方。

本实施例的方法在作用于投影图像的3步重建(数据加权，卷积滤波，反投影重建)之前增加了识别有效数据边界的步骤，即作为一种对投影数据的预处理，可用于和其他对反投影重建部分优化的算法联合使用，达到更理想的效果。

该方法只改变重建图像中背景，对重建物体的质量无任何影响。

在本发明实施例2中，举了一个应用实例，来具体阐述实施例1所述方法流程，实施例2具体包括以下步骤：

图4(a)为海螺图，是CT采集的海螺的3D投影数据中的一张切片。运用Ostu的方法求得该图像的最优二值化阈值后，对图4(a)进行二值化，二值化结果如图4(b)所示。由图4(b)，容易得到投影图中的有效数据边界，如图4(c)所示。

图5所示为对比应用优化前后重建方法所得到的重建结果。图5(a)、5(b)、5(c)为原始算法的重建结果。5(d)、5(e)、5(f)为采用本发明实施例的重建结果。其中，(a)、(d)为3D重建图像在y-z平面上的一个切片；(b)、(e)为3D重建图像在x-y平面上的一个切片；(c)、(f)为3D重建图像。并且，实验表明，由原始算法重建出图3(c)所需时间为运用本发明实施例中方法重建出图3(f)所需时间的5.2倍。图5中0.6cm是标尺，代表图上该长度的线段在现实中为0.6cm。

现有的与本发明实施例提供的重建方法的重建质量对比如图6所示。图6(a)、6(b)分别为图5(b)、(e)的灰度直方图。由于背景数据的灰度值小于90，所以在图6(a)、6(b)中，为了滤除背景数据，本发明实施例选择灰度值统计范围为90-255。图6(c)为图5(b)、5(e)白线处的灰度分布曲线，图中未加点的线代表图5(b)白线处的灰度分布曲线，加点的线代表图5(e)白线处的灰度分布曲线。图6(a)、6(b)表明，两种重建算法在整幅图像中灰度值范围为90-255的灰度直方图完全相同。这说明两种方法对目标物体的重建质量完全一样。图6(c)表明，两种重建算法的重建质量仅在背景处有所不同。而在白线(图5(b)、5(e)的白线)上的目标物体处，两种重建算法的重建质量完全一样。

为了找到速度提高倍数与物体尺寸的关系，本发明实施例重建了图7所示的一个椎体模型。图7(a)为椎体模型的3D重建效果图，图7(b)为x-y平面上的一张重建切片。

图8所示为本发明实施例对重建速度提高的效果。图8(a)描述图7中每张切片的重建速度，即重建速度与目标物体宽度的关系。图8(b)描述速度提高倍数与目标物体和探测器宽度之比的关系。图8(a)、8(b)中有几个奇点：A、B、F。A、F点的重建时间比旁边的稍慢是因为它们是第一张重建切片，在重建它们之前计算机需要将数据和代码读入内存。C点表示当物体在探测器上的的投影宽度和探测器尺寸接近时，本发明实施例的重建速度跟原始算法的重建速度没有太大区别。但是，在D点，可以看出当宽度比接近1.3时，重建速度可以提高到原来的2倍。这即说明，当物体投影宽度相当于探测器宽度的80％时，重建速度即可提高为原来的两倍。并且，当目标物体投影宽度为探测器宽度一半时，重建时间仅需原来的1/4。图8(c)为曲线拟合的结果。曲线方程为：y＝0.8811*x²+0.2427*x-0.1420。在该表达式中，二次项系数接近1并且远大于一次项和常数项的系数。一次项的系数不为0可能是由于在数据加权和卷及滤波时，对所有投影数据做了运算(并不是只对有效数据运算)。常数项系数不为0可能是与内存调度和MATLAB的编译方式有关。总的来说，这个拟合结果和上面理论分析的结果十分吻合。可以说明，速度提高倍数和目标物体投影宽度与探测器宽度之比为二次方的关系。

同时，由于该椎体模型每张切片的内部结构不尽相同，图8还说明了本发明实施例的提速效果与物体内部结构的复杂程度无关。

图9为本发明实施例提供的CT图像重建装置的结构图，该装置具体包括：

边界识别模块91，用于识别图像的有效数据边界，获得有效数据边界点行(列)号的最小值与最大值；

数据加权模块92，用于对图像进行数据加权，将扇形束校正为对应的平行束；

滤波模块93，用于通过卷积滤波消除图像中的伪影；

重建模块94，用于利用行(列)号的最小值与最大值进行反投影重建。

该装置的边界识别模块获取并保存图像中行(列)号的最小值与最大值，可以让计算机只处理行(列)号的最小值与最大值之间的数据(即有效数据)，从而提高计算机的重建速度。

其中，边界识别模块91具体可以包括：

读取单元911，用于读取所述图像数据；

数字处理单元912，用于采用数字图像处理的方法处理所述图像，寻找所述图像的像素值突变点；

返回单元913，用于返回所述像素值突变点行(列)号的最小值与最大值。

在本实施例中，该边界识别模块中的数字处理单元为二值化处理单元，采用二值化方式进行数字图像处理，实际上，该单元也可以采用边缘检测处理等方法进行数字图像处理，当进行二值化处理时，该单元还具体可以包括：

阈值获取器9121，用于获取所述图像的二值化最优阈值；

像素值设置器9122，用于遍历所述图像中的每个像点，如果当前像点的像素值大于或等于所述二值化最优阈值，则将此像点的像素值设置为1，反之则设置为0；

突变点获取器9123，用于将行(列)号相同像点的像素值累加，组成一个一维矩阵，将矩阵中第一个和最后一个值为非0的点设置为像素值突变点。

利用以上装置，通过识别有效数据的边界，之后仅对边界内的有效数据进行重建，在保证重建图像质量不变的前提下，使重建速度得到了提高。通过实验验证，重建速度具体提高的倍数为探测器尺寸与物体在探测器上投影最宽处尺寸之比的二次方。

总之，以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (8)

1、一种CT图像重建的方法，其特征在于，包括：

识别图像的有效数据边界，获得所述有效数据边界点行(列)号的最小值与最大值，对所述图像进行数据加权，将扇形束校正为对应的平行束，并通过卷积滤波消除所述图像中的伪影，利用所述行(列)号的最小值与最大值进行反投影重建。

2、根据权利要求1所述的方法，其特征在于，所述识别有效数据边界具体包括：

读取图像数据，采用数字图像处理的方法处理所述图像，寻找所述图像的像素值突变点，返回所述像素值突变点行(列)号的最小值与最大值。

3、根据权利要求2所述的方法，其特征在于，所述数字图像处理包括：

二值化处理或边缘检测处理。

4、根据权利要求3所述的方法，其特征在于，所述二值化处理的步骤具体包括：

获取所述图像的二值化最优阈值，遍历所述图像中的每个像点，如果当前像点的像素值大于或等于所述二值化最优阈值，则将此像点的像素值设置为1，反之则设置为0；

将行(列)号相同像点的像素值累加，组成一个一维矩阵，则矩阵中第一个和最后一个值为非0的点为像素值突变点。

5、根据权利要求2至4中任意一项所述的方法，其特征在于，所述对图像进行数据加权具体为：

对所述图像上的所有点进行数据加权或对所述有效数据边界点进行数据加权。

6、一种CT图像重建的装置，其特征在于，包括：

边界识别模块，用于识别图像的有效数据边界，获得所述有效数据边界点行(列)号的最小值与最大值；

数据加权模块，用于对所述图像进行数据加权，将扇形束校正为对应的平行束；

滤波模块，用于通过卷积滤波消除所述图像中的伪影；

重建模块，用于利用所述行(列)号的最小值与最大值进行反投影重建。

7、根据权利要求6所述的装置，其特征在于，所述边界识别模块具体包括：

读取单元，用于读取所述图像数据；

数字处理单元，用于采用数字图像处理的方法处理所述图像，寻找所述图像的像素值突变点；

返回单元，用于返回所述像素值突变点行(列)号的最小值与最大值。

8、根据权利要求7所述的装置，其特征在于，所述数字处理单元为二值化处理单元，具体包括：

阈值获取器，用于获取所述图像的二值化最优阈值；

像素值设置器，用于遍历所述图像中的每个像点，如果当前像点的像素值大于或等于所述二值化最优阈值，则将此像点的像素值设置为1，反之则设置为0；

突变点获取器，用于将行(列)号相同像点的像素值累加，组成一个一维矩阵，将矩阵中第一个和最后一个值为非0的点设置为像素值突变点。

Images