本发明涉及图像重建领域,具体涉及一种相对平行直线CT感兴趣区域图像重建方法。
背景技术:
最近提出的一种新的基于X射线源和探测器沿不同的方向平行移动(PTCT)的CT系统结构,并且已经证实它在低成本CT扫描仪有巨大的潜力。然而,为了优化PTCT系统,相关的图像重建应该被重点研究。
之前针对PTCT系统下FBP算法,主要处理完备且无截断的数据下的图像重建[3]。然而,在PTCT系统中探测器经常只能覆盖了物体的一部分,这导致数据被截断,进一步导致精确图像重建复杂,甚至不能进行图像重建。
技术实现要素:
有鉴于此,本发明的目的在于提供一种相对平行直线CT感兴趣区域图像重建方法。
本发明的目的是通过以下技术方案来实现的,
一种相对平行直线CT感兴趣区域图像重建方法,包括以下步骤:S1.获取投影数据;S2.对投影数据进行微分反投影;S3.对反投影获得的数据进行希尔伯特逆变换,获得重建图像。
进一步,所述步骤S2.反投影步骤是在线型PI线上产生一个中间希尔伯特图像函数其中
λ示射线源到原点的向量于y轴的夹角,h是原点到X射线源轨迹的距离,SDD为射线源轨迹到探测器轨迹的距离,ψ为射线源轨迹和x轴的夹角,λb和λe是X射线源轨迹开始和结束位置的角,p指投影数据。
进一步,在步骤S2中,希尔伯特逆变换得到的一次线性扫描所获得的图像为:
或
其中L>l≥max(xe|,|xb|),k(L,l,x)表示如下
进一步,所述l=max(|xb|,|xe|)+(2~3pixels),L=(1.1~1.3)max(|xe|,|xb|)。
进一步,从希尔伯特图像中获得的真实图像为:
ε是一个极小值,取值范围为(10-3,10-2)。
进一步,进行多次线性扫描,则在多次线性扫描模式下重建图像为:
或
进一步,还包括步骤S3.利用权重函数来削减多次直线扫描模型产生的冗余信息,所述权重函数为:
进一步,为了避免的不连续性,则
是一个光滑的正函数,
由于采用以上技术方案,本发明具有以下优点:
本发明提出的不仅能在截断不是很严重的情况下重建近似完整的物体图像,而且能够用PTCT系统收集的截断数据中无截断伪影精确重建感兴趣区域图像。
附图说明
为了使本发明的目的、技术方案和优点更加清楚,下面将结合附图对本发明作进一步的详细描述,其中:
图1为PTCT系统一般的直线扫描模型;
图2为PTCT系统数据获取模型;
图3为在平行直线扫描下给定方向通过终点关于不同扫描轨迹的投影数据;
图4为探测器截断感兴趣区域图像重建;
图5为用于重建的Sheep-Slogan图和笑脸图;
如图6分别使用公式(8)的FBP算法,公式(14),公式(26)的MP-BPF算法,和公式(27)的MZ-BPF算法对Shepp-Logan图进行无噪声扇形束非截断数据的图像重建,第一行到第三行分别代表1次,2次和3次扫描模型;
图7为2次扫描模型有权重函数和无权重函数的重建结果对比;
图8为MP-BPF,MZ-BPF和FBP算法在不同的扫描模型下y=0直线上的图像剖面图的灰度值大小;
图9为公式(8)的FBP算法,公式(15),公式(27)的MP-BPF算法,和公式(28)的MZ-BPF算法对Shepp-Logan图进行无噪声扇形束非截断数据的图像重建;
图10为使用公式(8)的FBP算法,公式(14),公式(26)的MP-BPF算法,和公式(27)的MZ-BPF算法对Shepp-Logan图进行无噪声扇形束截断数据的图像重建,第一行到第三行分别代表1次,2次和3次扫描模型;
图11为感兴趣区域图像重建;
图12为图10中y=0方向的剖面图;
图13为有噪声扇形束截断数据的图像重建为;
图14为图13沿y=0方向的剖面图。
具体实施方式
下面将结合附图,对本发明的优选实施例进行详细的描述。
首先回顾PTCT模型。在PTCT中,目标不动,X射线源和探测器沿相反方向移动,如图1所示。现在,以目标中心作为笛卡尔坐标系统的原点。射线源轨迹表示如下
λ示射线源到原点的向量于y轴的夹角,h是原点到X射线源轨迹的距离,SDD为射线源轨迹到探测器轨迹的距离。ψ为射线源轨迹和x轴的夹角。在这种工作下,假设图像函数
现在,介绍了一个以射线源点为圆心的移动坐标系统。在这个系统中,两个单位向量如下
事实上,在PTCT系统中采用平板探测器。然而,本文只考虑扇形束扫描。探测器上的任何参数都用t表示。投影p(t,λ,ψ)是沿指定的X射线通过射线源轨迹上的点和探测器元素t点的线积分。因此在t点的物象函数投影p(t,λ,ψ)表示如下
表示从射线源点到物象点方向的单位向量,表示如下
根据上述方程,
λb和λe是X射线轨迹开始和结束位置的角。是重建点到X射线源扫描轨迹的距离。g(t)是斜坡滤波器。
对于多次相对线性扫描模式,方程修改如下
根据经典的BPF算法,目标函数由扫描轨迹的弦组成,物体能够通过这些弦获得。然而对于PTCT系统,因为所有的弦覆盖在相同的线段上,与一般扫描轨迹的经典BPF算法相似,可以认为这些特殊的弦为直线型PI线,并且带有线型PI线与目标函数的交叉点是支撑段。
BPF算法的反投影步骤在线型PI线上产生一个中间希尔伯特图像函数
为了简化上述方程,投影数据的微分G(t,λ,ψ)表示如下
由方程(5)得
将方程(11),(12)带入方程(10)
将方程(1),(8),(13)带入方程(9)
B.希尔伯特逆变换
f(x)为有限区间[xb,xe]上的光滑函数。是它的一维希尔伯特变换,希尔伯特变换对如下
pv是积分的柯西主值。因为目标在有限区间[xb,xe],真实图像能够使用有限希尔伯特变换方程获得。在本发明中,采用了两个有限希尔伯特逆变换公式,分别为
其中L>l≥max(|xe|,|xb|),k(L,l,x)表示如下
其中f(ψ,λb,λe,x)是f(x)在一次线性扫描所获得的大致图像,理论上不是精确图像。在方程(18)中,对真实图像敏感的参数是L和l。实际中,可以在一个范围下选择一个适当的值l=max(|xb|,|xe|)+(2~3pixels),L=(1.1~1.3)max(|xe|,|xb|)
事实上,能够使用这些方程从希尔伯特图像中获得真实图像。然而,这些方程有一个奇点,为了避免积分的奇点,方程可以做如下修改
ε是一个极小值,取值范围为(10-3,10-2)。
基于获得的结果,在多次线性扫描模式下能够通过BPF算法被重建
上面两公式分别命名为MP-BPF算法和MZ-BPF算法。
现在,假定有N次直线扫描轨迹。为了方便,仅仅考虑第i次和第j次直线扫描轨迹,如图3所示。定义如下。是一个权重,目的是为了计算在多次线性扫描模型中获得的数据集的冗余信息。显然,在不同的扫描轨迹下当射线通过终点扇形束投影会为目标函数的重建提供重复的信息。
p(ti,λi,ψi)=p(tj,λj,ψj)s.t 1≤i,j≤N (23)
权重用来处理冗余信息。
是线与所有直线扫描路径的交点数,对于扫描轨迹上通过交叉点的固定直线为Φ(ψi,λi),因此,权重函数定义如下:
选择作为权重函数。然而,实际上,为了避免的不连续性,采取以下方法。
是一个光滑的正函数,本文表示为
最后,分别用MP-BPF算法和MZ-BPF算法实现图像重建。
总的来说,获得了从一系列扇形束平行直线扫描数据中重建图像的一般BPF算法。
在上述讨论中,仅仅考虑了这样一种情形,就是整个物体都能被探测器覆盖(无截断投影)。然而,很容易遇到这种情况。在PTCT系统中,目标物体不能被视野完全覆盖。(投影是截断的),如图5所示。
直线型PI线由参数组(ψi,λb,λe)表示,正如公式(26)和(27),要重建感兴趣区域图像仅仅需要知道在[xb,xe]和[-L,L]上的反投影图像。虽然不能从一次扫描模型中获得反投影图像,但是可以通过增加扫描线段的方法来获得完全投影实现完整图像的重建。例如二次和三次扫描模型,如果扩展这个结论,使用BPF算法精确重建图像必要的和充足的条件为:主函数上的任何点要被表示至少需要180度的扫描角。如果感兴趣区域被扫描线段完全照亮,可以使用本发明提出的BPF算法实现感兴趣区域图像重建,并且没有截断伪影。
对于扇形束PTCT系统提出的MP-BPF和MZ-BPF算法通过一个改进的Shepp-Logan图和一个笑脸图评估发现了它的有效性,如图6所示。笑脸图中被红色长方形覆盖的区域就是感兴趣区域,它被选择用来测试MP-BPF和MZ-BPF算法在截断投影下与FBP算法的比较的表现。两个图像的像素都是256×256它们覆盖的区域大小为1×1mm2。对于其它的1次,2次,3次扫描的参数如表格1所示。通过1000mm长的探测器阵列分别在1次,2次,3次扫描模式下产生非截断的PTCT数据。截断的数据由350mm长的探测器阵列获得。也通过高斯噪声加到无噪声的数据中来产生带噪声的数据。为了证明明显地低对比区域,选测无噪声数据的最大值的0.37%作为高斯噪声的标准差。将本文提出的BPF算法和FBP算法对非截断和截断数据分别做完全图像重建和感兴趣区域图像重建做对比。
表1仿真参数
如图6分别使用公式(8)的FBP算法,公式(14),公式(26)的MP-BPF算法,和公式(27)的MZ-BPF算法对Shepp-Logan图进行无噪声扇形束非截断数据的图像重建。第一行到第三行分别代表1次,2次和3次扫描模型。为了证实本文提出的权重函数的有效性,图7中分别给出了2次扫描模型有权重函数和无权重函数的重建结果对比。图8给出了MP-BPF,MZ-BPF和FBP算法在不同的扫描模型下y=0直线上的图像剖面图的灰度值大小。为了做鲜明的对比,真是图书的剖面图灰度值也呈现在图8中。为了进一步评估重建图像的质量,表2中列出了各算法的均方误差。
表2 PTCT系统使用FBP和BPF算法的均方误差
从以上分析可知,明显的可以通过2次和三次扫描获得精确的重建图像,3次扫描的图像质量优于2次扫描的图像质量。造成这种结果的主要原因是通过3次扫描收集的投影数据和圆形扫描模型一样绕物体扫描了360度,这里权重函数是一个简单的常数1/2。从另一个角度看,2次扫描模型被认为是一个典型的短扫描,因此权重函数很难确定。但是,本文提出的权重函数能够在一定程度上处理这些冗余信息,如图7所示。另外,在三次扫描模型下,对于无截断投影数据BPF算法的质量要差于FBP算法,因为使用BPF算法会对图像进行面远重置,这会降低图像质量。最后,结果证明MP-BPF算法和MZ-BPF算法在抑制图像噪声方面有很好的表现。
已经对笑脸图片应用了本文提出的BPF算法来重建图像的感兴趣区域。在本部分,需要观察完全近似图像获取,如图10所示,分别使用公式(8)的FBP算法,公式(14),公式(26)的MP-BPF算法,和公式(27)的MZ-BPF算法对Shepp-Logan图进行无噪声扇形束截断数据的图像重建。第一行到第三行分别代表1次,2次和3次扫描模型。
如图所示,可以看出使用BPF算法可以获得无截断伪影的重建图像,相比FBP算法就有截断伪影。另外,如果投影不是很严重截断的话,BPF算法能够重建超出感兴趣区域范围的完整的粗糙图像。
为了进一步研究本文提出的BPF算法在感兴趣区域图像重建方面的表现。如图11分别展示了从图10中获得感兴趣区域图像。为了比较感兴趣区域图像重建的质量,也在图12中提供了图11中y=0位置的剖面图。
可以看到,在图11中用红色长方形标注的截断伪影是使用FBP算法进行感兴趣区域重建时产生的。相比,使用BPF算法却可以没有截断伪影的精确重建感兴趣区域。在1次扫描模型下,因为投影的丢失,三种算法都不能重建感兴趣区域图像。在2次扫描模型下冗余信息是一个难题,另外使用BPF算法在3次扫描模式下可以获得比2次扫描模式下更完整的感兴趣区域重建图像。
为了进一步探索本文提出的PTCT系统的BPF算法对抑制图像噪声中的能力。如图13使用公式(8)的FBP算法,公式(14),公式(26)的MP-BPF算法,和公式(27)的MZ-BPF算法对Shepp-Logan图进行有噪声扇形束截断数据的感兴趣区域图像重建。图14是y=0方向的剖面图。这些图形证明本文提出的算法在抑制图像噪声方面有很好的表现。
最后说明的是,以上优选实施例仅用以说明本发明的技术方案而非限制,尽管通过上述优选实施例已经对本发明进行了详细的描述,但本领域技术人员应当理解,可以在形式上和细节上对其作出各种各样的改变,而不偏离本发明权利要求书所限定的范围。