CN100524358C - 一种改进的锥形束ct环形伪影的消除方法 - Google Patents

Abstract

本发明公开一种改进的锥形束CT重建图像环形伪影的消除方法，该方法将读入的CT重建图像先采用坐标变换的方法将环形伪影变成直线，再用与伪影平行的一维滤波器为掩模，对所述的中间图像进行逐步、逐行搜索，获取样本空间，同时计算出每一样本空间的像素值的方差，并根据该方差是否小于预定阈值来识别伪影，将样本空间内的每一像素点的像素值用该点及其纵轴方向上下两伪影之间像素点的平均像素值来替换从而消除环形伪影，最后将消除伪影后的中间图像变回到原坐标系下得到校正图像。

Description

一种改进的锥形束CT环形伪影的消除方法

技术领域

本发明涉及一种医学CT图像的处理方法，具体涉及在锥形束CT重建图像中环形伪影的消除方法。

背景技术

计算机断层成像技术，即CT(Computed Tomography)在20世纪临床医学领域得到了广泛的应用。CT自诞生以来，其扫描方式从二维的单一探测器平行束，发展到多探测器扇形束旋转扫描，三维CT则从最早的单排螺旋CT发展为多排探测器扫描，同时给出8-32层数据的多层面螺旋CT(Multislice spiral/helical CT，MSCT)。近几年使用面阵探测器的锥形束CT(CBCT)也逐渐进入了实用阶段，它采用锥形束扫描，获得各视角下的二维投影数据，重建后直接得到三维体积。相对于其它的几何结构，锥形束CT系统具有空间分辨率高、数据采集时间短、射线利用效率高等显著的优点。

按扫描轨迹来划分，CBCT的重建算法可以分为两大类：一类是Feldkamp、Dewis、Kress提出的锥顶轨迹为圆的重建方法(简称FDK算法)，另一类是Tuy、Smith、Grangeat提出的扫描轨迹为双圆或其它非平面轨迹的三维Radon变换方法。如果按算法的数学类型来划分，CBCT的重建算法可以分为分析算法和迭代算法。迭代算法比分析算法更占用计算机的系统资源，不满足实际工作中对重建速度必须是高速的要求，因此通常采用分析的算法。分析的算法又可以分为精确重建算法和非精确的重建算法。锥顶轨迹为圆时，不满足Tuy-Smith条件，其算法是一种非精确重建方法。但该方法具有扫描过程易于实现、计算相对简单的优点，因此得到了比较好的应用。基于单圆锥顶轨迹的FDK算法在radon空间不能得到完整的数据，当锥角比较小的时候缺失的数据区域很小，因此可以把这部分数据当作零或做插值来补齐这部分数据，然后再进行重建，可以得到比较好的近似结果。但是轨道半径变小，锥角增大，Radon空间中缺失的数据增多，这样圆轨道锥形束重建的近似效果就会下降，因此重建后的图像会出现明显的伪影。CBCT重建后的图像伪影主要有条形伪影(streak artifact)和环形伪影(ringartifact)，其中环形伪影是指以重建后图像等中心为中心的圆环形的伪影，环形伪影的主要是由X线探测器各个探测通道的响应不一致引起的。

CBCT重建后图像环形伪影的消除，主要包括环形伪影位置的检测和伪影信号强度的修正，即首先从重建后图像中定位出环形伪影的位置，然后根据伪影信号强度与其周围像素值大小对比，校正对应的伪影像素值。

环形伪影与其他类型的图像伪影相对比，是一种很有规律性的图像伪影。环形伪影的圆圈是以图像中心作为圆心的，这个性质使得在环形伪影定位的时候，可以先把图像从直角坐标变换到极坐标，再以极点为原点，极角为横轴，极径为纵轴构建一正交坐标系，然后将原图像的各个象素点的像素值分别变换到所述的正交坐标系中，原来的环形伪影在新的坐标系下就是一条直线。如此变换后，使用自定义的滤波器对所得到的直线伪影进行检测和替换处理，最后再将消除伪影的中间图象变换回原直角坐标下即可。GE医疗系统环球技术有限公司的专利申请(公开号：CN1720861A)公开了一种基于上述原理的环形伪影的消除方法。该方法定义的是纵轴方向的滤波器，且用滤波器下7个象素点中两头的两个点象素值和的2.5倍减去中心5个点象素值和的“差”来与系统设定上下限进行比较，进而判断滤波器中心象素是否有伪影，因此当伪影线的宽度大于或等于7个象素点的宽度时，所计算出的“差”就不落在系统设定的上下限内，必然会造成误判。此外，专利公开号为CN1720861A的专利申请方案仅用滤波器下7个点的象素值的平均值来校正伪影点的象素，那么当伪影线的宽度大于或等于7个象素点的宽度时，显然相当于没作替换。

发明内容

本发明的目的是提供一种改进的锥形束CT重建图像环形伪影的消除方法，该方法不受环形伪影宽度大小的约束，且不会破坏原有的图像信息。

本发明实现上述目的的技术解决方案是：

一种改进的锥形束CT重建图像环形伪影的消除方法，该方法包括以下步骤：

1)读入锥形束CT重建后图像，先以该图像所在直角坐标系的原点为极坐标系的极点，直角坐标系的横轴为极坐标系的极轴，直角坐标系中任一象素点到原点的距离为极径建立一极坐标系，再将所读入锥形束CT重建后图像转入建立的极坐标系中；

2)基于步骤1的极坐标系，以极角为横轴，极径为纵轴构建一正交坐标系，并按图像的每一象素点的极坐标，将原图像的各个象素点的像素值分别变换到所述的正交坐标系中，得到中间图像，在该中间图像中原直角坐标系中的环形伪影便转变为平行于横轴的直线；

3)在所述的正交坐标系中构造一横轴方向的一维滤波器，以该滤波器为掩模，对中间图象进行逐步、逐行搜索，获取样本空间，同时计算出每一样本空间的像素值的方差，并进行如下修正：

当所计算的样本空间的像素值的方差小于预定阈值时，将该样本空间内的每一象素点的像素值用该点及其纵轴方向上下两伪影之间象素点的平均像素值来替换，当所计算的样本空间的像素值的方差大于预定阈值时，该点的像素值不变；

其中以该滤波器为掩模对中间图像进行逐步搜索的步长小于或等于所述一维滤波器的长度；

4)修正完毕后，先以所述正交坐标系的原点为极点，横坐标为极角，纵坐标为极径将修正后的中间图像变回至极坐标系中，再以极坐标系的极点为原点，极角为横轴，任一象素点的极径为该点至原点的距离建立一直角坐标系，然后将修正后的中间图像转入所建立的直角坐标系中即可得到消除环形伪影后的校正图像。

为了便于计数，进一步提高处理速度和精度，本发明方法还可以包括下述步骤：在读入锥形束CT重建后图像后，建立极坐标系将图像的每一象素点的灰度值转换成浮点类型像素值；在得到消除环形伪影后的校正图像后，将图像中每一象素点的浮点类型像素值变换成像素值为0～255的灰度图像。

由于所述的环形伪影往往是非连续的，而待处理的环形伪影的直径大小又是不等的，因此所述的一维滤波器的矩阵大小不仅涉及到数据的处理量，还直接影响处理的精度。本发明人的实践证明，所述的一维滤波器的矩阵尺寸以1×3～9为宜。

对医学图像稍有研究者都知道，正常图像的灰度值在任一方向都是非线性的，且离散程度较大，而环形伪影变换为平行于横轴的直线后，水平方向上的它灰度值虽不是线性的，但离散程度小，本发明方法正是基于这一基本规律作出的技术方案。本发明上述方法中所述的样本空间的像素值的方差正是反映所述的离散程度的统计量，如果该统计量很小就说明其掩模的空间为环形伪影。因此，所述的预定阈值的大小直接影响伪影识别的准确性，本发明人在大量统计结果显示，所述的预定阈值为0.001～0.006可满足要求。

本发明的另一要点是伪影的消除。由于本发明的消除方法是将伪影点的像素值采用自身及其纵轴方向上下两伪影之间象素点的平均像素值来替换的，因此具体计算时，取点过多就有可能触及相邻的伪影，过少不具代表性，甚至可能取的都是伪影的象素点(当伪影线较宽时)，视两条伪影的距离，一般取伪影点及其纵轴方向上下各1～5个点的像素值来求均较为合适。

本发明对现有技术进行了重大改进，采用与伪影平行的一维滤波器为掩模，对所述的中间图象进行逐步、逐行搜索，获取样本空间，同时计算出每一样本空间的像素值的方差，并根据该方差是否小于预定阈值来识别伪影，将样本空间内的每一象素点的像素值用该点及其纵轴方向上下两伪影之间象素点的平均像素值来替换从而达到消除环形伪影的目的，因此具有下列突出的优点和效果：

1、分辨率高，伪影识别准确，且不受伪影宽度的限制；

2、伪影点的像素值的替换不触及相邻的伪影，不会破坏原有的图像信息。

附图说明

图1本发明方法的程序流程图；

图2为包含环形伪影的锥形束重建图像；

图3为在直角坐标系下含有环形伪影的锥形束重建图像；

图4为由直角坐标图像变换成极坐标图像；

图5为由极坐标图像换成正交坐标图像；

图6为消除环形伪影过程示意图；

图7为消除环形伪影后的中间图象；

图8为由正交坐标中间图象变换成极坐标中间图象；

图9为由极坐标中间图象变换成直角坐标图象；

图10为含有环形伪影的老鼠的锥形束CT重建后图像；

图11为消除环形伪影的老鼠的锥形束CT重建后图像。

具体实施方式

下面结合附图和一套用FDK算法重建后老鼠的锥形束CT图像来详细描述本发明的方法。具体步骤如下所述：

步骤1：读入用经过了FDK算法重建之后的老鼠的锥形束CT图像，其大小为512×256，通过线性变换将灰度图像转换成浮点类型的图像，以适合做坐标变换。具体运算公式为：I-I_min/(I_max-I_min)，式中I为图像像素值，I_max和I_min分别为待处理图像灰度最大值和最小值，在本例中分别为233和0。

步骤2：进行坐标转换(参见图2～4)。以CT图像所在直角坐标系的原点o为极坐标系的极点o，直角坐标系的横轴x为极坐标系的极轴，直角坐标系中任一象素点P到原点o的距离r为极径建立一极坐标系，那么CT图像直角坐标系中任一象素点P的象素值P(x，y)即可表达为P(r，θ)。表达式P(r，θ)中的θ范围为0～2π，r范围为 $0~\sqrt{{512}^2+{256}^2}/2.$ 将CT图像由直角坐标按公式 $\mathrm{\θ}=\arctan \left(\frac{y}{x}\right),$ $r=\sqrt{x^2+y^2}$ 转换为极坐标后，再以极径r为纵坐标，以极角θ为横坐标建立一正交坐标系(x′，y′)，将极坐标系中的CT图象变换至正交坐标系(x′，y′)中，作成中间图像(见图5)。在正交坐标系(x′，y′)中，重建后老鼠的锥形束CT图像中的环形伪影(见图2)即变成直线伪影。

步骤3：基于所述的正交坐标系(x′，y′)，构造一矩阵大小为1×6的x′轴方向的一维滤波器(参见图6)，以该滤波器为掩模，对中间图像进行逐步、逐行搜索，获取样本空间。在本例中以该滤波器为掩模对中间图像进行逐步搜索的步长为1，根据计算精度的需要，本例中滤波器的移动步长可取长1～6的任一整数值。在逐步、逐行搜索获取样本空间的同时，逐一计算出每一样本空间的像素值的方差，具体计算方法如下所述：由数字图象处理知识可知，所述的滤波器模板底下任一象素点所对应的像素值可表达为Pk(x′，y′)，样本均值可按公式 $\stackrel{\‾}{P}=\frac{1}{W}\underset{k=1}{\overset{W}{\mathrm{\Σ}}}{P}_k(x\′,y\′)$ 计算，方差可按公式 $S^2=\frac{1}{W-1}\underset{k=1}{\overset{W}{\mathrm{\Σ}}}{(P_k(x\′,y\′)-\stackrel{\‾}{P})}^2$ 计算，其中P为样本均值，S²为方差，W为一维滤波器的长度，在本例中W＝6。当S²小于或等于给定阈值0.002时，即按公式 $P_k^{\′}(x\′,y\′)=\frac{1}{5}\underset{l=1}{\overset{2}{\mathrm{\Σ}}}{P}_k(x\′,y\′+l)$ 计算出平均像素值，并用该值替换原象素值进行修正(对模板底下所有像素值逐一进行修正)，其中

为伪影修正后的像素值，k＝1～6；当S²大于0.002时，则不进行修正。按上述方法进行逐点、逐行修正后，即在正交坐标系(x′，y′)下得到修正后的中间图象(见图7)。由公式 $P_k^{\′}(x\′,y\′)=\frac{1}{5}\underset{l=1}{\overset{2}{\mathrm{\Σ}}}{P}_k(x\′,y\′+l)$ 可见，本实施例中的平均像素值是伪影点及其纵轴方向上下各2个象素点的平均值，视两条伪影的距离可适当地增减。

步骤4：将正交坐标系(x′，y′)下的修正后的中间图象变换到直角坐标系(x，y)下。先以所述正交坐标系(x′，y′)的原点为极点，横坐标为极角，纵坐标为极径将修正后的中间图像变回至极坐标系中(见图8)，再以极坐标系的极点为原点，并按公式x＝rcosθ，y＝rsinθ将修正后的中间图象变换到直角坐标系(x，y)下，其中，r是正交坐标系(x′，y′)的纵坐标值，θ是正交坐标系(x′，y′)的横坐标值。由于上述处理过程中所采用的是浮点类型的图像，因此最后可使用公式 $(I-I_{\min }^{\′})\×255/(I_{\max }^{\′}-I_{\min }^{\′})$ 进行线性变换，把它转换成灰度图像，其中

和

分别为处理后浮点图像灰度最大值和最小值，在本例中分别为1和0(参见图9)。经过上述处理后，所读入的CT图像(见图10)中的环形伪影完全消除，丝毫见不到修饰的痕迹，获得了理想的技术效果(见图11)。

Claims (5)

1、一种改进的锥形束CT重建图像环形伪影的消除方法，该方法包括以下步骤：

1)读入锥形束CT重建后图像，先以该图像所在直角坐标系的原点为极坐标系的极点，直角坐标系的横轴为极坐标系的极轴，直角坐标系中任一象素点到原点的距离为极径建立一极坐标系，再将所读入锥形束CT重建后图像转入建立的极坐标系中；

2)基于步骤1的极坐标系，以极角为横轴，极径为纵轴构建一正交坐标系，并按图像的每一象素点的极坐标，将原图像的各个象素点的像素值分别变换到所述的正交坐标系中，得到中间图像，在该中间图像中原直角坐标系中的环形伪影便转变为平行于横轴的直线；

3)在所述的正交坐标系中构造一横轴方向的一维滤波器，以该滤波器为掩模，对中间图像进行逐步、逐行搜索，获取样本空间，同时计算出每一样本空间的像素值的方差，并进行如下修正：

当所计算的样本空间的像素值的方差小于或等于预定阈值时，将该样本空间内的每一象素点的像素值用该点及其纵轴方向上下两伪影之间象素点的平均像素值来替换，当所计算的样本空间的像素值的方差大于预定阈值时，该点的像素值不变；

其中以该滤波器为掩模对中间图像进行逐步搜索的步长小于或等于所述一维滤波器的长度；

4)修正完毕后，先以所述正交坐标系的原点为极点，横坐标为极角，纵坐标为极径将修正后的中间图像变回至极坐标系中，再以极坐标系的极点为原点，极角为横轴，任一象素点的极径为该点至原点的距离建立一直角坐标系，然后将修正后的中间图像转入所建立的直角坐标系中即可得到消除环形伪影后的校正图像。

2、根据权利要求1所述的一种改进的锥形束CT重建图像环形伪影的消除方法，其特征在于：

步骤1)在读入锥形束CT重建后图像后，建立极坐标系将图像的每一象素点的灰度值转换成浮点类型像素值；

步骤4)在得到消除环形伪影后的校正图像后，将图像中每一象素点的浮点类型像素值变换成像素值为0～255的灰度图像。

3、根据权利要求1或2所述的一种改进的锥形束CT重建图像环形伪影的消除方法，其特征在于所述的一维滤波器的矩阵尺寸为1×3～9。

4、根据权利要求1或2所述的一种改进的锥形束CT重建图像环形伪影的消除方法，其特征在于所述的预定阈值T为0.001～0.006。

5、根据权利要求1或2所述的一种改进的锥形束CT重建图像环形伪影的消除方法，其特征在于所述的平均像素值为样本空间内的伪影点及其纵轴方向上下各1～5个象素点的平均像素值。

Images