CN111487265A - 一种结合投影一致性的锥束ct硬化伪影校正方法 - Google Patents

Abstract

本发明提供一种结合投影一致性的锥束CT硬化伪影校正方法。该方法包括：根据多能投影对数化数据重建得到CT图像，对CT图像进行分割，并进行二值化，得到二值化图像模型；根据实际扫描时的锥束CT成像***参数建立重投影坐标系，将二值化图像模型重投影到平板探测器上，通过二值化图像模型和虚拟探测器之间的对应关系对X射线与二值化图像模型进行求交计算，得到射线穿越物体长度；根据射线穿越物体长度和多能投影对数化数据进行多项式拟合，得到硬化伪影校正模型；根据硬化伪影校正模型校正多能投影对数化数据，得到等效单能投影对数化数据；根据等效单能投影对数化数据重建得到校正后的CT图像。本发明无需任何前置条件，适用于多种类型探测器。

Description

一种结合投影一致性的锥束CT硬化伪影校正方法

技术领域

本发明涉及平板探测器射束硬化校正技术领域，尤其涉及一种结合投影一致性的锥束CT硬化伪影校正方法。

背景技术

在计算机断层成像(Computed Tomography，CT)技术中，对于能量为E₀、强度为I₀的单能射线，根据比尔定律，可知其穿过长度为l的物体后强度I为：I＝I₀exp(-∫u(x,y,z)dl)；其中，u(x,y,z)为被测物体某断面在E₀情况下线性衰减系数的分布函数。单能投影对数化数据p_m为：

量分布在(E_max,E_min)范围内，设其能谱分布为S(E)，且线性衰减系数为能量的函数，则由比尔定律可得

多能投影对数化数据p_p为：

由p_p的表达式可知，l越大，p_p就越大，但两者之间不是线性关系。设p_p＝g(l)，g为射线硬化模型，为讨论该函数的增减情况，对其求导可知，g'(l)是关于l的单调递减函数，即随l的增大p_p增加的幅度变小。由于射线束中低能光子较高能光子衰减得快，因此，在穿越不同厚度的物体时射线束的能量分布亦不同，射线束的平均能量对应的衰减系数u不再是常数，而是随穿越物体厚度的增加而减小。由于射线硬化效应，真实的对数化投影数据p_p与重建中假设的对数化投影数据p_m关于线性衰减系数u(x,y,z,E)具有不同的表达式，该数据与穿越长度不是线性关系，重建过程中用p_p代替p_m，从而在重建图像中造成伪影，这种因为多能X射线造成的伪影即为硬化伪影，表现为单一材质物体重建图像中心的灰度值较边缘小且呈现杯状。硬化伪影会严重降低CT图像质量，因此需要对其进行校正。

针对硬化伪影校正，目前已经有许多种校正方法，主要有多项式拟合法(黄魁东等.基于切片轮廓重投影的锥束CT射束硬化校正方法[J].仪器仪表学报.2008.029(009)：1873-1877)、蒙特卡罗校正方法、迭代校正法等。多项式拟合法的优点在于多项式系数容易计算，思路简单易于接受，校正效果在一定条件下较好。鉴于该方法的有效性和简便性，是目前最常用的校正方法。但多项式拟合的方法有放大噪声的作用，在校正射束硬化的同时也放大了图像中的噪声，且该方法对于单种材质的被检测物效果较好，但需要有与被检测物相同材质的体模，实际应用中缺乏灵活性。蒙特卡罗法也称统计模拟法，是一种以概率统计理论为指导的非常重要的数值计算方法，其核心是随机抽样方法。该方法应用的前提是必须已知被检测物体的材料。迭代法，是一种不断用变量的旧值递推新值的过程。基于迭代法的一个优点在于一些物理效应，如射束硬化可以融合到正投影的过程中，这样射束硬化就可以考虑到重建过程中去，以达到消除伪影的目的。但由于迭代重建算法计算量大、重建速度慢、可并行度低等原因，目前仍处于研究阶段，实际应用中并不广泛。

发明内容

为解决现有校正方法中因需要前置条件(例如需要与被检测物相同材质的体膜，需要提前知道被检测物体的材料)而灵活性较差以及计算量大的问题，本发明提供一种结合投影一致性的锥束CT硬化伪影校正方法。

本发明提供的结合投影一致性的锥束CT硬化伪影校正方法，包括：

步骤1：根据多能投影对数化数据重建得到CT图像，并对所述CT图像进行图像分割，对分割后得到的轮廓图像进行二值化，得到物体实体区域的二值化图像模型；

步骤2：根据实际扫描时的锥束CT成像***参数建立重投影坐标系，将所述二值化图像模型重投影到平板探测器上，通过二值化图像模型和虚拟探测器之间的对应关系对X射线与二值化图像模型进行求交计算，得到射线穿越物体长度；

步骤3：根据所述射线穿越物体长度和所述多能投影对数化数据进行多项式拟合，得到硬化伪影校正模型；

步骤4：根据所述硬化伪影校正模型校正所述多能投影对数化数据，得到等效单能投影对数化数据；

步骤5：根据所述等效单能投影对数化数据重建得到校正后的CT图像。

进一步地，步骤2具体为：

首先确定被重建点(x,y,z)在中心平面的投影点C，然后从C点分别作辅助线CA、CD垂直于X轴与中心射线SO，得到几何关系：

FA＝xtanβ，CF＝y-xtanβ

CE＝CF×cosβ＝-xsinβ+ycosβ

又由于：

a＝R×tanγ

故射线在虚拟探测器上的水平位置a(x,y,β)为：

进而得到射线在虚拟探测器上的竖直位置b(x,y,z,β)为：

其中，O表示旋转中心位置，(x,y,z)表示被重建点坐标，S表示X射线源位置，so′是锥形X射线的中心射线，SK表示经过被重建点的一条射线，k′代表K点在平板探测器中心层面的投影，k′的被重建点为M，U为SM长度，β为中心射线与Y轴的夹角，k为射线SK的锥角，γ为sk′与中心射线so′的夹角，E、F分别为中心射线SO与CD、CA的交点；A、D均为X轴上的点；R表示X射线源到物体中心的距离。

进一步地，步骤4包括：

步骤4.1：给定实验物质和X射线对应所述实验物质的透射厚度x，对所述实验物质进行扫描得到X射线透射强度I_p对应的射线和数据y：

x＝[x₁ x₂…x_n]^T

y＝ln(I₀/I_P)＝[y₁ y₂…y_n]^T

步骤4.2：根据所述透射厚度x与所述射线和数据y，采用曲线y＝ax^b对所述射线和数据y进行拟合，利用最小二乘法拟合估计得到曲线中a、b的参数a′和b′：

步骤4.3：根据拟合估计得到的a′和b′，计算等效透射厚度x_eq：

步骤4.4：根据所述等效透射厚度x_eq，得到等效射束和比尔定律的拟合方程；

步骤4.5：将多能投影对数化数据带入步骤4.2得到的拟合曲线

中，得到与所述多能投影对数化数据相对应的穿透厚度x1，再将所述穿透厚度x1带入步骤4.4得到的等效射束和比尔定律的拟合方程中，得到与所述多能投影对数化数据相对应的等效单能投影对数化数据。

本发明的有益效果：

本发明提供的结合投影一致性的锥束CT硬化伪影校正方法，通过利用射线贯穿物体长度与多能投影值之间的关系曲线，再对该曲线进行拟合的相关关系进行射束硬化校正，首先对原始投影的重建图像进行预处理，并进行图像分割得到图像分割信息；然后根据实际扫描时的锥束CT成像***参数建立重投影坐标系,通过射线与二值化图像模型重投影到平板探测器上的求交计算，得到射线穿越物体长度与投影灰度之间的对应关系；最后拟合射束硬化曲线并校正射束硬化伪影，并利用数据一致性条件进行优化。本发明无需前置条件(例如无需与被检测物相同材质的体膜，也无需提前知道被检测物体的材料)，所需硬化校正参数均由对投影图像进行计算所得，因此可适用于多种探测器，并且相对于迭代重建算法，本发明的计算量较小，可简单有效地消除锥束CT的射束硬化伪影。

附图说明

图1为本发明实施例提供的一种结合投影一致性的锥束CT硬化伪影校正方法的流程示意图；

图2为本发明实施例提供的扫描结构示意图；

图3为本发明实施例提供的圆周扫描几何关系示意图；

图4为本发明实施例提供的射束硬化伪影校正前的图像；

图5为本发明实施例提供的利用本发明方法进行射束硬化伪影校正后的图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

如图1所示，本发明实施例提供一种结合投影一致性的锥束CT硬化伪影校正方法，该方法包括以下步骤：

S101：根据多能投影对数化数据重建得到CT图像，并对所述CT图像进行图像分割，对分割后得到的轮廓图像进行二值化，得到物体实体区域的二值化图像模型；

S102：根据实际扫描时的锥束CT成像***参数建立重投影坐标系，将所述二值化图像模型重投影到平板探测器上，通过二值化图像模型和虚拟探测器之间的对应关系对X射线与二值化图像模型进行求交计算，得到射线穿越物体长度；

S103：根据所述射线穿越物体长度和所述多能投影对数化数据进行多项式拟合，得到硬化伪影校正模型；

S104：根据所述硬化伪影校正模型校正所述多能投影对数化数据，得到等效单能投影对数化数据；

S105：根据所述等效单能投影对数化数据重建得到校正后的CT图像。

本发明实施例提供的结合投影一致性的锥束CT硬化伪影校正方法，不需要任何的前置条件(例如无需与被检测物相同材质的体膜，也无需提前知道被检测物体的材料)，所需硬化校正参数均由对投影图像进行计算所得，因此可适用于多种类型探测器，并且相对于迭代重建算法，本发明的计算量较小，可简单有效地消除锥束CT的射束硬化伪影。

实施例2

在上述实施例1的基础上，本发明实施例还提供一种结合投影一致性的锥束CT硬化伪影校正方法，包括以下步骤：

S201：根据多能投影对数化数据重建得到CT图像，并对所述CT图像进行图像分割，对分割后得到的轮廓图像进行二值化，得到物体实体区域的二值化图像模型；

具体地，图像分割时，对重建得到的CT图像分割地越精细越好，至少应将被测物体的大致轮廓分割出来。重建后的CT图像生成了一个立体数据，为了获得每一个体数据所代表的信息而进行的二值化，二值化之后可以清楚地获得每一个体数据的信息，例如0代表空气，1代表物体。实体区域就是物体所在的区域。

S202：根据实际扫描时的锥束CT成像***参数建立重投影坐标系，将所述二值化图像模型重投影到平板探测器上，通过二值化图像模型和虚拟探测器之间的对应关系对X射线与二值化图像模型进行求交计算，得到射线穿越物体长度；

具体地，将所述二值化图像模型重投影到平板探测器上，得到二值化图像模型和虚拟探测器之间的对应关系，包括：

圆轨迹FDK重建算法扫描结构如图2所示，其中O表示旋转中心位置，(x,y,z)表示被重建点坐标，为了方便表示与计算，在O点位置引入了虚拟探测器，被扫描物***于扫描中心位置。X射线源位于S，so′是锥形X射线的中心射线，SK表示经过被重建点的一条射线，k′代表K点在平板探测器中心层面(Z＝0)的投影，其被重建点为M，SM长度为U，中心射线与Y轴成β角，k为射线SK的锥角，γ为sk′与中心射线so′的夹角。

圆周扫描几何关系如图3所示，图中C为坐标(x,y,z)的被重建点在中心平面的投影点，SO为中心射线，从C点分别作辅助线CA、CD垂直于X轴与SO，可以得到如下几何关系：

FA＝xtanβ，CF＝y-xtanβ

CE＝CF×cosβ＝-xsinβ+ycosβ

又由于：

a＝R×tanγ

故射线在虚拟探测器上的水平位置a(x,y,β)为：

进而得到射线在虚拟探测器上的竖直位置b(x,y,z,β)为：

具体地，根据以上公式推导可以获知二值化图像模型和虚拟探测器之间的对应关系，从而得到单个虚拟探测器元素所对应的二值化图像模型元素的个数，可以表示为X射线穿透的二值化图像模型的个数，即射线穿越物体长度。

S203：根据所述射线穿越物体长度和所述多能投影对数化数据进行多项式拟合，得到硬化伪影校正模型；

S204：根据所述硬化伪影校正模型校正所述多能投影对数化数据，得到等效单能投影对数化数据；

具体地，在校正过程中，首先将通过实验测定一组透射厚度x的射束和数据y，拟合后得到X射线射束和数据y与透射厚度x之间的关系曲线，即X射线射束和拟合曲线。然后，利用X射线射束和拟合曲线拟合出等效透射厚度，从而得到等效射束和比尔定律的拟合方程。最后再进行重建，就可有效消除X射线射束硬化的影响。作为一种可实施方式，本步骤包括以下子步骤：

S2041：给定实验物质和X射线对应所述实验物质的透射厚度x，对所述实验物质进行扫描得到X射线透射强度I_p对应的射线和数据y：

x＝[x₁ x₂…x_n]^T

y＝ln(I₀/I_P)＝[y₁ y₂…y_n]^T

其中，n表示X射线穿透物体部位不为0的个数；

S2042：根据所述透射厚度x与所述射线和数据y，采用曲线y＝ax^b对所述射线和数据y进行拟合，利用最小二乘法拟合估计得到曲线中a、b的参数a′、b′：

具体地，采用曲线y＝ax^b对射线和数据y进行拟合的方法是：两边取对数得

lgy＝lga+blgx

则a′，b′可用最小二乘法拟合估计：

S2043：根据拟合估计得到的a′和b′，计算等效透射厚度x_eq：

具体地，X射线射束和拟合曲线可用拟合方程：

来拟合。对应地，多色射束和的透射厚度x欲修正为单色射束和的等效透射厚度，则需保证X射线等效为单色射线的衰减系数与射束和值，由比尔定律可知，透射厚度x应表为等效透射厚度x_eq，即等效射束和比尔定律应表为y＝ux_eq，从而可得到x_eq。

S2044：根据所述等效透射厚度x_eq，得到等效射束和比尔定律的拟合方程；

S2045：将多能投影对数化数据带入步骤S2042得到的拟合曲线

中，得到与所述多能投影对数化数据相对应的穿透厚度x1，再将所述穿透厚度x1带入步骤S2046得到的等效射束和比尔定律的拟合方程中，得到与所述多能投影对数化数据相对应的等效单能投影对数化数据。

S205：根据所述等效单能投影对数化数据重建得到校正后的CT图像。

为了验证本发明提供的锥束CT硬化伪影校正方法的有效性，本发明还提供了图4和图5：图4为未利用本发明方法进行校正的重建图像；图5为利用本发明方法进行校正后的重建图像。如图4所示，可以看出图像上存在严重的杯状伪影，射束硬化现象严重。而如图5所示，经过本发明方法校正后，本发明对射束硬化造成的杯状伪影能达到很好的去除效果，相比于传统校正方法对杯状伪影校正得更加彻底，硬化伪影校正效果良好。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (3)

1.一种结合投影一致性的锥束CT硬化伪影校正方法，其特征在于，包括：

步骤1：根据多能投影对数化数据重建得到CT图像，并对所述CT图像进行图像分割，对分割后得到的轮廓图像进行二值化，得到物体实体区域的二值化图像模型；

步骤2：根据实际扫描时的锥束CT成像***参数建立重投影坐标系，将所述二值化图像模型重投影到平板探测器上，通过二值化图像模型和虚拟探测器之间的对应关系对X射线与二值化图像模型进行求交计算，得到射线穿越物体长度；

步骤3：根据所述射线穿越物体长度和所述多能投影对数化数据进行多项式拟合，得到硬化伪影校正模型；

步骤4：根据所述硬化伪影校正模型校正所述多能投影对数化数据，得到等效单能投影对数化数据；

步骤5：根据所述等效单能投影对数化数据重建得到校正后的CT图像。

2.根据权利要求1所述的方法，其特征在于，步骤2具体为：

首先确定被重建点(x,y,z)在中心平面的投影点C，然后从C点分别作辅助线CA、CD垂直于X轴与中心射线SO，得到几何关系：

FA＝x tanβ，CF＝y-x tanβ

CE＝CF×cosβ＝-x sinβ+y cosβ

又由于：

a＝R×tanγ

故射线在虚拟探测器上的水平位置a(x,y,β)为：

进而得到射线在虚拟探测器上的竖直位置b(x,y,z,β)为：

其中，O表示旋转中心位置，(x,y,z)表示被重建点坐标，S表示X射线源位置，so′是锥形X射线的中心射线，SK表示经过被重建点的一条射线，k′代表K点在平板探测器中心层面的投影，k′的被重建点为M，U为SM长度，β为中心射线与Y轴的夹角，k为射线SK的锥角，γ为sk′与中心射线so′的夹角，E、F分别为中心射线SO与CD、CA的交点；A、D均为X轴上的点；R表示X射线源到物体中心的距离。

3.根据权利要求1所述的方法，其特征在于，步骤4包括：

步骤4.1：给定实验物质和X射线对应所述实验物质的透射厚度x，对所述实验物质进行扫描得到X射线透射强度I_p对应的射线和数据y：

x＝[x₁ x₂ … x_n]^T

y＝ln(I₀/I_P)＝[y₁ y₂ … y_n]^T

其中，n表示X射线穿透物体部位不为0的个数；

步骤4.2：根据所述透射厚度x与所述射线和数据y，采用曲线y＝ax^b对所述射线和数据y进行拟合，利用最小二乘法拟合估计得到曲线中a、b的参数a′、b′：

步骤4.3：根据拟合估计得到的a′和b′，计算等效透射厚度x_eq：

步骤4.4：根据所述等效透射厚度x_eq，得到等效射束和比尔定律的拟合方程；

步骤4.5：将多能投影对数化数据带入步骤4.2得到的拟合曲线

中，得到与所述多能投影对数化数据相对应的穿透厚度x1，再将所述穿透厚度x1带入步骤4.4得到的等效射束和比尔定律的拟合方程中，得到与所述多能投影对数化数据相对应的等效单能投影对数化数据。

Images