基于光子计数探测器的CT系统在蒙特卡罗模拟平台的实现
2018-09-10孙永刚孔慧华张海娇
孙永刚 孔慧华 张海娇
摘 要:为了解决利用Matlab,C++等软件对投影图像、重建算法的模拟无法真实反映CT系统的成像过程的问题,在研究模拟理论的基础上,首先利用Geant4和Gate数据包搭建蒙特卡罗模拟平台;然后通过光子计数探测器模拟多能谱CT成像,同时分析在不同标准下的阈值划分对重建图像的影响,利用光子计数探测器同时扫描出多个不同能段下的投影图像,并对比出各个能段下的投影效果;其次,根据不同能段的划分,选取每个能段的360度下的投影数据,利用FBP算法重建投影切片图像;最后,对所得到的切片图像进行材料区分,得到硫和铝的衰减系数曲线。实验结果表明,模拟方法对光子计数探测器的成像过程做出了模拟,得到了较好的模拟图像,并给出了相关阈值的划分依据,可为其他CT成像研究提供参考。
关键词:图像处理;多能谱CT;能谱划分;Geant4;Gate
中图分类号:TP391.4 文献标志码:A
文章编号:1008-1534(2018)05-0335-06
可见光的差别在于波长和频率,在各个能量下波长和频率是不同的,因此可以利用这一现象(即光谱探测原理),设计具有X射线能量分辨能力的多能谱CT投影系统[1]。陈平等[2]基于能谱滤波分离获取多个能谱通道下的投影。杨晓飞[3]设置2个射线源得到2个能谱投影图像,并对光子计数型探测器设定阈值获取多能谱投影。能谱分离一是通过加装有滤波片的探测器,识别不同能量下的光子,以此实现多能谱的获取;二是设置多个探测源,可以利用多个入射能量获得不同能量的投影。当前多数情况下主要利用双源双探测器实现对2个不同能谱的分辨。但上述办法效率不高,投影过程相对繁琐;三是先设定好需要的能量阈值,然后探测器可以区分X射线的光子脉冲,对选定的不同能谱区域内的光子数进行累计计数,从而获取多个能谱通道下的投影。但光子计数探测器的模拟成本高昂,步骤复杂,因此搭建计算机模拟平台成为理想的选择。上述3种能谱获取方式如图1所示。
笔者通过搭建蒙特卡罗平台来模拟光子计数探测器的成像过程,通过Geant4数据包的数据算法程序计算探测器接收到的光子个数,然后通过Gate平台的模拟程序,设计出光子计数探测器,在进行实验后,得到了每个能谱阈值之间的投影测量数据;最后利用投影测量数据用FBP算法重建出切片图像[4],以此分析不同能谱通道的划分对图像的影响。
1 光子计数探测器成像模型
光子计数探测器能对X射线中的不同能量的光子分别记录,可获得多个不同能谱通道下的投影,根据兰博比尔(Lambert-Beers)定理,一个阈值S下的光子强度为
2 蒙特卡罗模拟平台
Geant4可对X射线中不同粒子的运动过程进行仿真模拟,具有开放且可编辑的源代码,可根据具体的实验需求构建出CT仿真系统。Gate是用于对CT成像的模拟软件,可以设置对应的探测器的参数,比如探测器像素、材料及厚度。同时设置好模体的材料及形状,最后取得对模体扫描后的投影图像,进而完美地实现对CT系统的模拟仿真[5]。
2.1 蒙特卡罗模拟平台的搭建
平台主要利用Geant4和Gate数据包进行模拟实验,2种数据包需要在Linux系统上运行。为了便于仿真模拟,笔者在Windows系统环境下,通过Red Hat公司开发的Cygwin虚拟镜像软件,配置运行环境,搭建出模拟平台的外部框架。安装VC++2010进行环境变量的设置,并在安装完成后运行实例程序进行检测。
至此平台的搭建工作基本完成,在进行模拟实验的时候,通过Gate调用Geant4数据包,实现对CT成像过程的模拟。在Gate中,可以分别设计CT系统的点源、模体及探测器,同时可以根据具体实验需求来设定所需的入射能量、物理过程及数据输出的种类。
2.2 平台对X射线散射过程的模拟
在一个CT系统里,首先要设置相应的X射线源,在蒙特卡罗模拟平台里,可以通过点源文件设置入射光子强度、入射角度以及点源数量。为了提高模拟效率,可以设置一个锥束X射线源,射线源发出的X射线到达模体后,会产生一系列物理现象,例如光电效应、韧致辐射、康普顿散射等。
在模拟康普顿散射及韧致辐射现象时,Geant4数据包具有多种散射模型,例如所包含的G4-KleinNishinaCompton,G4-KleinNishinaMode模型。第1种模型的模拟机制不同,其模拟速度较快,第2种模型中考虑了粒子的各种物理现象,包括原子壳效应、光电效应,对比之前的方法更为客观,但模拟速度也随之变慢,因此笔者采用第2种散射模型。所模拟的散射过程符合“克来茵-仁科”截面公式:
平台中的韧致辐射模型可以对低能时的入射截面进行计算,对应的X射线成像过程也可获得准确的模拟。每次仿真模拟时,首先要确定每个光子的入射能量,入射角度及光子个数,然后对光子的所有物理过程进行模拟,记录下光子的运动信息。当入射的光子产生物理反应时,根据康普顿、韧致辐射截面公式,可以计算出光子的随机概率分布,同时根据一定规则对入射光子进行筛选,从而得到了到达探测器的光子个数N及光子能量E,最后由此生成投影的測量数据。
2.3 扫描模体的设置
笔者选取一个包含了4个不同半径且不同材料的子圆柱体的有机圆柱体,4个子圆柱体的材料分别为金属物质铝和铁,非金属物质碳和硫,外部的母圆柱体材料选取人体骨骼有机化合物。其横截面切片图如图2所示。模体中的具体参数如表1所示。
2.4 光子计数探测器的设计
当前光子计数探测器选取的晶体材料大多为硅(Si),砷化镓(GaAs),碲化镉(CdTe)和碲锌镉(CZT),不同探测材料的能量分辨率和探测范围也不尽相同。考虑材料的性能及成本,设置厚度1.5 mm,大小1.0 mm×1.0 mm的碲锌镉(CZT)材料为探测器的晶体单元。这里构建的面阵探测器包含256×256个探测单元,设置圆柱体载物台,系统旋转中心轴为载物台中心轴。射线源高度与载物台中心轴处于同一平面,这样就构成一个锥束CT成像系统,见图3。
3 模拟实验
3.1 基于光子计数型的多能谱CT模拟
平台通过根据CT系统的各部分来设置对应的模块文件,以此完成模拟实验,模块文件主要有X射线源模块Source.mac,扫描模体模块Phantom.mac以及探测器主体模块CTScanner.mac。
在Source.mac模块文件中,设置好入射点源位置、入射强度及粒子种类;在Phantom.mac模块中设置扫描模体的几何形状、模体材料和模体坐标;在CTScanner.mac中设置探测器材料及像素大小。在平台上模拟出CT成像过程,得到对应的投影数据文件[7]。
设置好X射线源、扫描模体、光子计数探测器及输出文件后,通过SpectrumGUI软件进行能谱数值模拟,获取能谱为0~140 keV、电压在125 kVP时的离散能谱数据,然后加入到source.mac文件中,即可得到模拟所需要的X射线源的能谱分布数据文件。首先模拟截取能谱为10~20 keV,25~35 keV,60~90 keV,95~125 keV下的投影过程,每个通道下进行360度的扫描,并得到所有能谱通道的投影数据;然后选取角度在137°时的投影数据,并生成投影图像;最后利用Matlab軟件,由FBP算法重建出投影图上第69层的切片图像。能谱图像及对应结果如图4和图5所示。
分析重建好的切片图像,通过对比可以看出不同能谱通道下的投影数据也有不同,根据徐品等[8]的研究理论,计数值在最低时的阈值对应较高的入射电压,所以在高电压下光子计数探测器的光子计数值在较高能量下会相对较少,探测器在高能量下接收到的光子数减少,使得部分物质的重建亮度降低。为进一步验证实验结论,需要做出每个能谱通道下所对应的衰减系数曲线图像,可以选取部分材料区域,以便对衰减系数与理论数据作出对比。其中深色是实验数据,浅色是理论数据(该数据选取平台模拟不添加生成噪声的物理现象,但受其他因素诸如探测器的探测效率及模拟硬件的限制,噪声无法完全避免,只能选取噪声最低时的效果做对比),结果如图6所示。
可以看出不同的扫描有不同的衰减系数曲线,在重建出的切片图像中,图像颜色越浅,该区域的衰减曲线越接近0。图像的噪声越小,曲线波动也就越小,而密度较高的元素变化越大。低能谱时的光子强度比高能谱小,轮廓分明。对比硫和铝区域的衰减系数曲线图可以得到,在25~35 keV之间,曲线更为平缓,与理论数值的拟合较好,因此图像质量更好。但分别在25~35 keV,60~90 keV和95~125 keV之间时,当光子数减少时,图像对比度在60~90 keV,95~125 keV之间比在25~35 keV之间要小,衰减曲线较之前波动大,与理论数值拟合较差,物质成分的区分度降低。
3.2 能谱划分对投影的影响
在得到多个能谱下的投影后,可以看出,在不同能段下,探测器对光子的吸收效率不同。而在相同的能谱范围里,无法保证探测器在每一段能谱下也能接受相同数量的光子数,这样会对所生成的投影图像有一定影响。
为此,设置不同的阈值,但要使探测器能接收数量大致相同的光子数;同样再设置范围相同的阈值做对比。阈值设定范围在16 keV左右,投影、重建结果所对应的衰减系数如图7所示[9]。
在探测器接收的光子数接近时,投影结果在高能和低能时比较接近,重建出来的切片图像中,金属与非金属物质区别较大,不同材料的投影对比度也有较大的差异,但衰减系数的曲线更加平稳。在选取相同阈值范围时,低能谱下的投影对比度有较大差异;在高能时,探测器能接收到的光子数减少,图像物质识别度降低。
因此,只有合理设定能谱划分阈值,使探测器能在每个能谱通道接受数量相同的光子数,才能获得更好的投影效果[10]。
4 结 论
通过研究多能光子计数探测器的成像过程,首先验证了在蒙特卡罗模拟平台可以准确模拟CT成像;然后对相应的物理过程、能谱阈值的选择及成像效果有了更深层次的研究,得到的投影图像真实可靠,更加容易识别出不同材料的物质;最后根据实验得到的图像,分析物质的衰减变化,充分验证了实验结果的准确性。
利用计算机平台模拟CT系统的成像过程,丰富了相关的研究方法,进一步简化了实验过程,降低了实验所需要的人力、物力等成本,可为其他的CT成像研究提供模拟手段及研究基础。
参考文献/References:
[1] 庄天戈.CT原理和算法[M].上海:上海交通大学出版社,1992.
[2] 陈平,郭蓉,潘晋孝, 等. 基于Geant4的能谱滤波分离虚拟平台研究[J]. 中国测试,2016,42(6):65-69.
CHEN Ping,GUO Rong,PAN Jinxiao,et al. Study on virtual platform based on Geant4 energy spectrum filtering separation[J].China Measurement & Test, 2016,42(6):65-69.
[3] 杨晓飞. 基于GATE平台的双能CT成像系统模拟设计[D].沈阳:东北大学, 2014.
YANG Xiaofei. The Simulation Design of Dual Energy CT Imaging System Based on GATE Platform[D]. Shenyang: Northeastern University, 2014.
[4] 郭蓉,潘晋孝,陈平,等. X射线谱的Geant4模拟及其滤波效应研究[J]. 核电子学与探测技术,2015,35(8):767-769.
GUO Rong, PAN Jinxiao, CHEN Ping, et al. Geant4 simulation of X-ray spectrum and the study of filtering effect[J]. Nuclear Electronics & Detection Technology,2015,35(8):767-769.
[5] 廖玉婷,王永波,许可欣,等. 基于光子计数探测器的X线CT蒙特卡罗仿真[J]. 中国医学物理学杂志 2016,33(2):122-127.
LIAO Yuting, WANG Yongbo, XU Kexin, et al. X-ray computed tomography Monte Carlo simulation based on photon countingdetector[J]. Chinese Journal of Medical Physics,2016,33(2):122-127.
[6] Geant4 Collaboration.Geant4 User's Guide for Application Developers(Version:geant4 10.0)[R].[S.l.]:CERN,1996:159-163.
[7] BAUMER C, MARTENS G, MENSER B, et al. Testing an energy-dispersive counting-mode detector with hard X-rays from a synchrotron source[J]. IEEE Transactions on Nuclear Science, 2008,55(3): 1785-1790.
[8] 徐品,陳奭,袁刚,等.基于光子计数探测器的能谱CT的研究[J].CT理论与应用研究, 2015,24(6):822-823.
XU Pin, CHEN Shi, YUAN Gang, et al. Study of spectral CT based on a photon counting detector[J]. CT Theory and Applications, 2015,24(6):822-823.
[9] 徐亚东, 介万奇, 查刚强, 等. CdZnTe平面探测器对低能X/γ射线的光谱响应[J]. 光学学报, 2009, 29(11): 3072-3077.
XU Yadong, JIE Wanqi, ZHA Gangqiang, et al. A study on the low energy X/γ-ray spectral resonse of CdZnTe planar detectors[J]. Acta Optic Sinica, 2009, 29(11): 3072-3077.
[10] 陈平, 潘晋孝, 刘宾. 连续能谱X-CT投影仿真算法[J]. 无损检测, 2009, 31(2): 102-104.
CHEN Ping, PAN Jinxiao, LIU Bin. Simulation arithmetic of X-CT projection based on consecutive spectrum[J]. Nondestructive Testing, 2009, 31(2): 102-104.