基于图形处理器的并行加速双能计算机层析成像重建方法 
下载： 720次
1 引言
对于计算机层析成像(CT)重建领域,采用图形处理器(GPU)并行加速已经成为重要的方法,具体的加速方法既是本领域研究人员研究的核心技术,也是提高CT产品实用性的关键,直接决定了产品的市场竞争力。现有文献中已有多种针对单能CT重建的GPU加速方法。与单能CT重建不同,双能CT重建并不是前者的简单重复,而是具有其特殊性。
现有文献中双能CT重建算法大致可以分为3类:迭代重建算法[1-2]、后处理重建算法[3-5]和预处理重建算法[6-9]。由于迭代重建方法计算速度慢,不适合在对处理速度要求较高的安检CT设备中使用,目前安检双能CT设备主要采用后两种算法——后处理重建算法和预处理重建算法,这2种重建方法各有特点。后处理重建方法能够获得较高的形状精度,适合展示给用户观察,但由于不能从根本上消除射束硬化效应,重建的被检物材料精度难以保证,不适用于自动材料识别。预处理重建算法理论上可以得到不受能谱影响的被检物的等效原子序数和电子密度信息,重建的材料精度较高,适合用于自动材料识别,但是图像形状精度损失严重,不适合用户观察。双能CT要想同时获得上述后处理重建和预处理重建图像的各自优势,就需要重建4种图像,如何准确和快速地重建图像是双能CT能否用于实际应用的关键。
目前在双能CT的快速准确重建实现上还存在一些制约因素,现有文献中有多种双能投影分解方法,如多项式拟合法[6]、等值线法[7]、高次曲面函数拟合法[9-10]、迭代求解法[11-13]、优化的分解方法[14]、基于投影匹配的投影分解算法[15]等,为了在实际应用中快速准确地进行投影分解,采用双能投影分解查找表方式。但采用这种方式能够保证投影分解准确性的前提是输入查找表的每对双能投影数据穿过物体的路径必须相同。近年来,采用GPU进行并行加速成为CT重建加速的研究热点,CT的各种重建算法基本都实现了GPU加速[16-22]。但是对于双能CT重建加速,现有文献没有给出有针对性的有效解决方案,仅仅是通过简单重复执行单能CT重建加速方法来实现双能CT重建加速[20],重建2种图像的时间翻倍,重建4种图像的时间将至少是单能CT重建时间的4倍,这将严重影响实际应用。
本文提出了一种基于GPU的并行加速双能CT重建方法,该方法采用上述后处理的线衰减系数重建和预处理的基材料(基效应)重建算法同时重建4种图像,获得了各自重建图像的优势,有效解决了双能CT重建的精度和速度问题,在实际应用中取得良好的效果。
2 GPU并行加速双能CT重建方法
2.1 双能CT重建流程
双能CT扫描一般采用两种方式获得双能采样——真双能和伪双能。真双能采用两种能量不同的射线进行两次扫描获得采样数据;伪双能采用单一能量的射线源,利用夹心探测器获得两个能量的采样,夹心探测器通过在两个探测器之间安装滤波片实现双能探测[23]。真双能采样通过扫描机械控制精度误差,存在高低能采样数据视角不一致的问题。双能CT需要同时采用后处理的线衰减系数重建和预处理的基材料(基效应)重建生成4种重建图像。如果只是采用后处理的线衰减系数重建方法,对两个能量分别采用单能CT重建算法逐个进行重建,即使高低能采样数据视角不一致,也只是使重建的高低能线衰减系数图像整体转了一定的角度,不会导致重建图像数值错误。但对于预处理的基材料(基效应)重建,如果存在高低能采样数据视角不一致的问题则会导致双能投影分解错误,从而导致后续重建结果出现错误。为了保证重建的准确性,本文提出的双能CT重建流程先要进行高低能采样数据的视角配准过程,以保证后续双能投影分解的准确性,重建流程如
双能CT重建流程首先进行双能采样数据的视角配准,然后再计算原始的双能投影数据,后处理的线衰减系数重建方法可直接对原始的双能投影数据进行投影加权、滤波和反投影重建;而预处理重建方法则需要对原始的双能投影数据进行投影分解,获得基材料(或基效应)的投影数据,然后再对分解(基材料或基效应)投影进行投影加权、滤波和反投影重建,最后再根据基材料(或基效应)的重建结果计算等效原子序数图像和电子密度图像。其中双能采样数据的视角配准的方法采用文献[
16]的基于双能正弦图求解最大相关系数的配准方法,双能采样数据视角配准流程如
图 2. 双能采样数据视角配准的流程图
Fig. 2. Flow chart of view angle registration of dual-energy sample data
视角配准流程先对采集的双能采样数据的左部分(或右部分)进行二值化,并将二值化后的图像进行横向投影获得二值图的一维投影信号,再计算双能一维投影信号的互相关系数数组以确定相关系数数组中的最大值所在位置,从而得到双能采样数据的相对视角差。
2.2 算法并行性分析和GPU加速方法设计
通过对2.1节所述双能CT重建流程进行并行性分析,设计基于GPU的并行算法。基于参考文献[ 23],双能采样数据的视角配准是通过计算双能采样数据f0(x,y)与f1(x,y)的一维投影Pf0(y)与Pf1(y)的最大归一化互相关系数来实现,归一化互相关系数ρ(Pf0,Pf1)的计算公式为
式中:n为一维投影Pf0(y)与Pf1(y)的尺寸;A、B、C、D均为常数,可分别表示为
式中:
假设采样数据的宽度(即探测器单元数量)为w,高度(即采样视角数量)为h。为加快处理速度,算法实现时仅对有效宽度为v(v<w)、高度为h的采样数据进行二值化,并计算二值化图像的一维投影Pf0(y)与Pf1(y)。在GPU加速实现时,可以通过一个内核函数按照视角并行同时实现高低能采样数据的二值化和高低能二值化图像的一维投影计算,具体过程如
配准后高低能采样数据的投影计算、投影分解及投影加权都可以按照射线并行合并在一个内核函数中以实现加速,从而得到原始高低能投影的加权结果和分解后的基材料(或基效应)投影的加权结果。合并在同一内核函数中实现上述3步计算的优点在于节省存取中间结果的时间,加快处理速度。其中双能投影分解基于投影分解查找表实现,预先将离线生成的投影分解查找表存储到二维纹理存储器中以提高读取速度,该二维纹理采用非归一化的浮点型拾取坐标,将高低能投影值作为查找时输入的横纵坐标,并利用纹理存储器硬件实现的线性模式滤波功能对读取的浮点型返回值进行插值,从而实现快速高精度的双能投影分解。
通过频域点积实现卷积滤波步骤。采用CUDA提供的快速傅里叶变换CUFFT函数库,在初始化时先通过快速傅里叶变换将滤波核数据变换到频域,并将其存储到一维纹理存储器中,以便在频域运算时加快读取速度。投影数据进行FFT变换时采用复数形式,利用复数的实部和虚部同时进行双能加权后投影的频域变换及分解加权后投影的频域变换,并将FFT后的复数值存储到复数类型的纹理存储器中以便计算频域点积时加快读取速度。
反投影重建步骤可以按照像素并行在同一内核函数中进行4种图像的反投影重建,即一个线程完成4种类型图像(即高、低能线衰减系数图像和2种分解图像)的同一像素位置的反投影重建计算,参见
反投影重建实现时也可以通过将单个重建像素点的全部投影视角数平均分配到多个线程,从而在实现像素并行的同时实现视角维度的并行计算。假设设计用M个线程并行处理单个像素点的全部视角叠加。每个线程仅完成X个均匀分布视角的反投影计算,X=N/M,N为扫描一圈采集的全部视角数量,M为用来并行处理视角的线程数量,且能够被N整除。内核函数计算完X个视角后再给全局存储器的相应像素位置赋值,以便节省反投影重建的访存时间。可以根据具体采用的GPU硬件环境的不同设置不同的M值, M等于1时,对重建的像素点不存在 “写冲突”,在将累加结果写入全局存储器时可以不使用原子加操作;当M大于1时,由于增加了按视角并行的线程,在将累加结果写入全局存储器时需要使用原子加操作。
重建流程的最后一步是根据之前重建的两种分解(基材料或基效应)图像计算等效原子序数图像和电子密度图像。具体实现时也可以按照像素并行在同一内核函数中进行计算,一个线程仅计算等效原子序数图像和电子密度图像上一个像素点。
完整的GPU并行加速的双能CT重建流程和并行模式如
图 3. 基于GPU并行加速的双能CT重建的流程和并行模式图
Fig. 3. Flow chart and parallel mode of parallel accelerated dual-energy CT reconstruction based on GPU
3 实验分析
为了评价该方法,在NVIDIA GeForce GTX TITAN Black显卡上采用宽和高分别为896(探测器单元数量)、720(采样视角数量)的双能CT采样数据进行重建尺寸为512 pixel×512 pixel的重建实验180次,不包括上述流程中的采样数据配准步骤,重建过程和其中最费时的反投影过程的平均耗时结果如
表 1. 不包括采样配准的重建过程总耗时及反投影步骤耗时
Table 1. Spending time of back-projection step and total reconstruction time without sampling registrationms
|
重建1幅图像(即高、低能线衰减系数图像)的180次重建总耗时均值为10.7759 ms,如果采用现有文献方法并通过简单重复执行来实现双能CT重建2幅图像(即高低能线衰减系数图像)的耗时至少为21 ms,而采用本文的方法重建2幅图像(即高低能线衰减系数图像)的重建总耗时均值为11.1696 ms,实验条件下总耗时的加速比为1.88。如果采用现有文献方法并通过简单重复执行实现双能CT重建4幅图像(即高、低能线衰减系数图像和2幅基材料分解的图像)的重建总耗时至少为42 ms,而采用本文方法重建这4幅图像的耗时均值为12.9669 ms,总耗时的加速比为3.24。采用本文方法重建4+2幅图像(即高、低能线衰减系数图像,2幅基材料分解的图像,以及由2幅基材料分解的图像再计算出等效原子序数图像和电子密度图像)的耗时均值为13.3755 ms,由于计算等效原子序数图像和电子密度图像的过程没有费时的反投影步骤,实验条件下计算过程耗时仅占重建4幅图像时总耗时的3%。
图 4. 双能CT重建图像。 (a)高能线衰减系数图像;(b)低能线衰减系数图像;(c)等效原子序数图像;(d)电子密度图像
Fig. 4. Reconstruction images of dual-energy CT. (a) High-energy linear attenuation coefficient image; (b) low-energy linear attenuation coefficient image; (c) effective atomic number image; (d) electron density image
在单能CT重建步骤中,反投影步骤是最耗时的,占重建总时间的80%以上,重建1幅图像的反投影耗时的180次均值为8.6519 ms,如果仅仅是通过现有文献中简单重复执行单能CT重建来实现双能CT反投影2幅图像,耗时至少为17 ms,重建4幅图像反投影耗时至少为34 ms,而采用本文方法重建2幅图像的反投影耗时的均值为8.9416 ms,反投影2幅图像的加速比为1.9,重建4幅图像的反投影耗时的均值为9.2883 ms,实验条件下反投影4幅图像的加速比为3.66。反投影耗时与重建过程总耗时的比值呈下降趋势。
4 结论
针对双能CT重建,提出了一种基于GPU的并行加速重建方法,与采用现有单能CT加速重建方法重复计算的两个能谱相比,该方法在高低能采样数据视角配准的前提下,同时重建4种图像,提高了重建的准确性,同时获得后处理重建图像的形状精度和预处理重建图像的材料精度优势,使用GPU并行加速方法实现算法全过程,进行快速配准;合并了原始投影计算、投影分解和投影加权的并行加速过程,在加权步骤避免了对原始双能投影和分解后投影权值的重复计算;在滤波步骤,通过并行FFT将投影数据和滤波核变换到频域中,通过频域点积实现滤波,充分利用复数的实部和虚部进行双能同步滤波;在反投影步骤同时重建4幅图像,节省了重复计算4幅重建图像的反投影地址的时间,同时也缩短了计算中间过程的数据存取时间。实验证明了本文方法的有效性,本文方法对双能CT成像的实际工程应用具有积极的意义。
[1] Fessler J A, Elbakri I A, Predrag S, et al. Maximum-likelihood dual-energy tomographic image reconstruction[J]. Proceedings of SPIE, 2002, 4684: 38-49.
[2] Sukovic P, Clinthorne N H. Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography[J]. IEEE Transactions on Medical Imaging, 2000, 19(11): 1075-1081.
[3] Christ G. Exact treatment of the dual-energy method in CT using polyenergetic X-ray spectra[J]. Physics in Medicine and Biology, 1984, 29(12): 1501-1510.
[4] Steenbeek J C M, van Kuijk C, Grashuis J L, et al. Selection of fat-equivalent materials in postprocessing dual-energy quantitative CT[J]. Medical Physics, 1992, 19(4): 1051-1056.
[5] Goodsitt M M. Beam hardening errors in post-processing dual energy quantitative computed tomography[J]. Medical Physics, 1995, 22(7): 1039-1047.
[6] Alvarez R E. MacOvski A. Energy-selective reconstructions in X-ray computerised tomography[J]. Physics in Medicine and Biology, 1976, 21(5): 733-744.
[7] Chuang K S, Huang H K. A fast dual-energy computational method using isotransmission lines and table lookup[J]. Medical Physics, 1987, 14(2): 186-192.
[8] Alvarez R E, Seibert J A, Thompson S K. Comparison of dual energy detector system performance[J]. Medical Physics, 2004, 31(3): 556-565.
[9] Cardinal H N, Fenster A. An accurate method for direct dual-energy calibration and decomposition[J]. Medical Physics, 1990, 17(3): 327-341.
[10] Cardinal H N, Fenster A. Analytic approximation of the log-signal and log-variance functions of X-ray imaging systems, with application to dual-energy imaging[J]. Medical Physics, 1991, 18(5): 867-879.
[11] Tang S J, Mou X Q, Luo T. An iteration algorithm in dual-energy X-ray imaging based on polychromatic physics model[J]. Proceedings of SPIE, 2006, 6142: 61422R.
[12] Zou Y, Silver M D. Analysis of fast kV-switching in dual energy CT using a pre-reconstruction decomposition technique[J]. Proceedings of SPIE, 2008, 6913: 691313.
[13] 李保磊, 张萍宇, 李斌, 等. X射线双能计算机层析成像投影分解的优化迭代方法[J]. 光学学报, 2017, 37(10): 1034001.
[14] Ying ZR, NaiduR, GuilbertK, et al. Dual energy volumetric X-ray tomographic sensor for luggage screening[C]∥2007 IEEE Sensors Applications Symposium, February 6-8, 2007. San Diego, CA, USA. IEEE, 2007: 1- 6.
[15] 李保磊, 张耀军. 基于投影匹配的X射线双能计算机层析成像投影分解算法[J]. 光学学报, 2011, 31(3): 0311002.
[16] Xu F, Mueller K. Accelerating popular tomographic reconstruction algorithms on commodity PC graphics hardware[J]. IEEE Transactions on Nuclear Science, 2005, 52(3): 654-663.
[17] Xu F, Mueller K. Real-time 3D computed tomographic reconstruction using commodity graphics hardware[J]. Physics in Medicine and Biology, 2007, 52(12): 3405-3419.
[19] 魏存峰. CT图像重建的GPU加速[D]. 北京: 中国科学院高能物理研究所, 2008.
Wei CF. Acceleration of CT image reconstruction[D]. Beijing: Institute of High Energy Physics Chinese Academy of Sciences, 2008.
[20] 毕文元. 双能量单层螺旋CT物品机的图像重建方法研究[D]. 北京: 清华大学, 2010.
Bi WY. Dual energy single slice helical CT reconstruction algorithms for luggage security inspection system[D]. Beijing: Tsinghua Univesity, 2010.
[21] 张慧滔, 于平, 胡修炎, 等. 利用GPU实现单层螺旋CT的三维图像重建[J]. 电子学报, 2011, 39(1): 76-81.
Zhang H T, Yu P, Hu X Y, et al. Single-slice helical CT three-dimensional image reconstruction using GPU[J]. Acta Electronica Sinica, 2011, 39(1): 76-81.
[22] 首都师范大学. CT图像重建的GPU加速方法: CN200810113846.0[P].2010. 12. 15.
Capital NormalUniversity. GPU acceleration method for CT image reconstruction: CN200810113846.0[P].2010. 12. 15.
[23] 李保磊, 张耀军. 针对X射线双能CT成像的正弦图快速配准方法[J]. 光学技术, 2011, 37(2): 198-202.
张萍宇. 基于图形处理器的并行加速双能计算机层析成像重建方法[J]. 激光与光电子学进展, 2020, 57(12): 121001. Pingyu Zhang. Parallel Accelerated Reconstruction Method for Dual-Energy Computed Tomography Based on Graphics Processing Unit[J]. Laser & Optoelectronics Progress, 2020, 57(12): 121001.
PDF全文
