基于折射角的稀疏投影角度相衬CT图像重建

司 凯,孙丰荣,宋尚玲,秦 峰

(1.山东大学信息科学与工程学院,济南250100;2.山东大学第二医院设备部,济南250033)

基于折射角的稀疏投影角度相衬CT图像重建

司 凯1,孙丰荣1,宋尚玲2,秦 峰1

(1.山东大学信息科学与工程学院,济南250100;2.山东大学第二医院设备部,济南250033)

已有基于X射线吸收衬度机制的计算机断层成像(CT)技术很难对由轻元素构成的弱吸收物质进行高质量成像。X射线相位衬度CT成像是对弱吸收物质具有超高分辨率的一种CT技术,但该技术成像时间长、所需X射线辐射剂量大,不利于临床推广应用,因此,研究稀疏投影角度条件下的X射线相位衬度CT图像重建问题,基于压缩感知图像重建理论,使用折射角信息减少X射线辐射剂量,提出一种X射线相位衬度CT图像重建算法。实验结果表明,与滤波反投影算法相比,该算法在稀疏投影角度下可以得到较高质量的重建图像,在实际数据实验中能获得较高的峰值信噪比和数值准确性。

X射线相衬;计算机断层成像;图像重建;压缩感知;折射角

1 概述

在1895年发现了X射线并拍摄了第一幅X光片后,随着计算机技术的日渐成熟,X射线计算机断层成像(Computerized Tomography,CT)技术应运而生,并在医学诊断、材料学、工业无损检测和安保等领域得到了广泛的应用。传统的X射线成像技术是基于物体各个部位对于硬X射线吸收量的不同来重建图像,但是由轻元素构成的弱吸收物质的吸收量很小,成像效果不理想,因而限制了它的进一步应

用。然而当X射线穿过物体时,它的强度和相移共同发生变化,在这个变化中本文使用δ代表X射线穿过物体发生的相移量,β则代表物体对X射线的吸收量。

研究发现,对于人体肺部组织等弱吸收物质,X射线穿过时相移量δ的变化等级可达吸收量β的千倍以上［1］,使用硬X射线穿过物体的相移量来成像可以达到高密度分辨率,因此基于相移量的成像技术无疑更加适应弱吸收物质。在此背景下,文献［2-3］分别于1995年和1996年在《Nature》上提出了硬X射线相位衬度(简称相衬)成像技术,开创了基于相移量δ成像技术的研究领域。

近年来,X射线相衬成像技术发展出4种主要的方法:干涉仪成像法［4］、衍射增强成像(Diffraction Enhanced Imaging,DEI)［5］法、类同轴成像法［6］和光栅干涉仪成像法［7］。本文方法使用折射角信息并针对稀疏投影角度进行图像重建,算法适应于衍射增强成像法和光栅干涉仪成像法,可以改善传统相衬成像技术的采集时间长、辐射剂量大等问题。

目前,关于X射线相衬CT图像重建技术的研究主要有2个方向:第1个方向是重建出相位量的梯度▽δ［8-9］;第2个方向是重建出相位量δ［10-11］。文献［12］分析将压缩感知(Compressed Sensing,CS)理论与重建相位量梯度结合的可行性,提出利用经典ART迭代方法求解该问题的算法,并利用实验验证了其算法的可靠性。文献［13］在求解相位量梯度的问题中,也得到了较好的实验结果。以上述研究为基础,本文使用折射角信息,基于恢复δ的重建策略,在CS图像重建理论的框架下,采用基于距离驱动(Distance Driven,DD)的正/反投影运算,并引入有序子集(Ordered Subset,OS)等思想,提出一种基于CS的相衬CT迭代图像重建算法。

2 CT图像重建中的压缩感知理论

在稀疏投影角度下的CT图像重建属于求解欠定线性系统的问题,即:

其中,A是m×n的系统矩阵,且m＜n;矩阵v∈Rn×1为待重建图像;m维列向量u为投影数据。结合CS理论,如果待重建图像v是S稀疏的,且测量矩阵A(CT图像重建中的系统矩阵)满足严格的等距同构特性(Restricted Isometry Property,RIP)［14］,即对于1≤S≤n,可以定义一个等距同构常数γS,满足:

可以将式(1)转化成下述有约束的l1范数最优化问题［15-16］:

其中,Ψv为图像v的稀疏化表示;Ψ为稀疏基矩阵。据此,当系统矩阵A满足RIP时可以通过求解式(3)得到欠定线性系统式(1)理论上的精确结果。在实际应用中,系统矩阵A很难满足RIP性质,并且存在信号只能达到近似稀疏或投影数据含有噪声等情况。虽然如此,鉴于求解有约束的l1范数最优化问题的鲁棒性,此领域的学者普遍将系统矩阵A近似的视作CS中的测量矩阵,并依然可以得到很好的重建效果［13,17］。

3 本文算法介绍

在本文算法中,为满足式(3)中对目标函数的稀疏性要求,本文将x的梯度做为目标函数。在迭代重建过程中,本文采用了距离驱动的正/反投影运算,引入了自适应性权重和OS思想并自适应地选择OS数量,使得算法的整体性能大大提高。在式(3)的最优化环节,可以采用任一种现有的寻优算法,本文中使用模拟退火(Simulated Annealing, SA)算法。

3.1 折射角预处理

本文采用重建δ项的策略,如图1所示,其中,l⊥代表射线传播的垂直方向;p-q为旋转坐标,将折射角定义为X射线穿过物体后发生偏折的角度θ,它与相位项的二维空间分布δ(x,y)的关系如下［2］:

图1 折射角光路图

在图1的光路图中,相位项δ满足以下等式［18］:

(1)将采集到的折射角信息θ(p,β)与符号函数sgn(p)做卷积的预处理操作,得到投影数(即式(1)中的y)。

(2)采用本文重建算法对式(3)中有约束的l1范数最优化问题进行求解。

3.2 启发式求解有约束l1范数最优化问题的方法

求解过程简要伪代码如下:

可以看出,迭代重建环节分为2个部分,粗略重建和最优化,分别对应本节启发式策略中步骤(1)、步骤(2)的内容。其中,粗略重建部分内部包含IterationⅡ和IterationⅢ2种迭代。在最优化部分,首先利用粗略重建后的图像跟初始图像计算偏差距离,然后使用SA算法得到最优解。

3.2.1 有序子集思想和自适应性权重

在IterationⅢ中,以往迭代类算法使用一个依据经验选择的松弛因子来促进算法收敛,但是松弛因子的选定比较繁琐,自适应性差。本文利用迭代次数和当前重建图像像素值提出一种自适应性的权重w(3.2节求解过程简要伪代码第20行)替代松弛因子,可以增强本文算法对于不同投影数据的适应能力。对于第m次内部外层迭代:

其中,t(r,c,d)为距离驱动的正投影运算中第(r,c)个像素对第d个检测器的贡献值。

3.2.2 距离驱动的正/反投影运算

在IterationⅢ中,现有的正/反投影运算主要分为3类,即基于像素驱动、基于射线驱动和基于距离驱动的方法。基于像素驱动的方法硬件实现简单,但是存在过高的运算复杂度以及在正投影运算中容易引入高频噪声等弊端。基于射线驱动的方法在反投影运算中会产生摩尔纹伪影,其运算方式会导致不连续的内存访问机制从而影响处理速度。而文献［20］提出的基于距离驱动的方法具有低算术复杂度、高度连续的内存访问机制和正/反投影运算中不易产生伪影等优点［21］。鉴于这3类运算方法的优缺点,本文采用基于距离驱动的正/反投影运算完成本文方法中的迭代环节,以期获得更高的算法通用性和高质量重建结果。

3.3 算法实现过程

基于上述分析以及旨在重建δ项的初衷,本文将折射角信息预处理后做为投影数据,使用OS思想

从时间向度看，任何国家和地区的社会管理模式与思想都不是一成不变的，都有一个生成与转换的过程。任何社会管理思想都深深依存于它所产生的社会环境，目的也是为了解决其当下社会面临的秩序、发展、社会经济文化等问题。从追求政治和行政独立，到以绩效为核心，再到更为重视社会公正，发展任务的转换决定了不同时期西方社会管理的特点。

结合距离驱动正/反投影运算的CT迭代图像重建技术得到粗略重建结果,然后通过SA算法得出最优解。综上,算法可分为两部分,折射角预处理、迭代重建,其中,迭代重建分为粗略迭代重建和最优化计算。算法描述如下:

(1)折射角预处理

将折射角信息利用3.1节中叙述的方法与符号函数进行卷积操作得到真实投影数据矩阵P,其角度数为N。

(2)迭代重建

假定待重建图像为矩阵V,算法总的迭代次数i=1,2,…,I。对于每一次总的迭代,由粗略迭代重建和最优化计算两部分组成:

1)粗略迭代重建

此过程内部包含2种迭代:一种是粗略迭代重建部分自身的迭代,表示为IterationⅡ(迭代次数m=1, 2,…,M);一种是每一次IterationⅡ下所有投影角度的迭代,表示为IterationⅢ(迭代次数n=1,2,…,N)。

每一次IterationⅢ,即在每一个投影角度下,根据3.2.1节中方法自适应地划分K个有序子集。对于第n次IterationⅢ,每个有序子集(k=1,2,…,K)都进行正向投影、作差并加权、反向投影以及图像修正4个过程,其中,算法初始值V(i=1,m=1,n=1,k=1)=0,具体步骤如下:

②作差并加权:将相应的真实投影数据P(i,k)与w×PDD_OS作差,得到差矩阵PDiff。

③反向投影:使用距离驱动的反向投影运算对PDiff进行反向投影得到偏差修正图像矩阵VDiff。

④图像修正:将本次IterationⅢ的初始图像与偏差修正图像作和,得到第n次IterationⅢ中第k个有序子集的重建图像,即V(i,m,n+1,k)=VDiff+V(i,m,n,k)。

2)最优化计算

最优化计算环节需要根据上述粗略迭代重建过程得出计算的范围,所以本环节可分为求解偏差距离和寻优两部分。具体如下:

①求解偏差距离:

令i=i+1,并将V(i)作为本次总体迭代的初始值,重复粗略迭代重建和最优化计算过程,直到i=I,程序终止,得到最终重建结果V=V(i=I)。

4 实验结果与分析

对于X射线相衬设备减小X射线辐射剂量最直接的方式就是稀疏投影角度扫描。本文分别使用仿真和真实数据的实验结果来验证稀疏投影角度下本文方法的性能,并将其与滤波反投影(Filtered Back-projection,FBP)算法结果进行比较。欠采样扫描是一种不满足奈奎斯特抽样定理的投影数据采集方式,本文对于X射线源从0°～180°旋转,主要使用3种欠采样投影数据采集方式:(1)180个角度模式,采集步长为1°;(2)90个角度模式,采集步长为2°; (3)36个角度模式,采集步长为5°。

本文使用以下3种通用图像质量评判标准来衡量本文方法的性能:

(1)归一化均方根误差(Normalized Root Mean Square Error,NRMSE):

其中,Xu,v和Yu,v分别为重建图像和参考图像在(u,v)时的像素值,图像的大小为(U,V)。NRMSE的值越小,说明重建图像有越高的数值精确性。

(2)峰值信噪比(Peak Signal to Noise Ratio, PSNR):

其中,Xmax为重建图像中像素值的最大值。PNSR值越大,说明重建图像失真越小。

(3)统一质量指标(Universal Quality Index, UQI),用于衡量2个图像接近程度的指标,其值越高,重建图像越接近参考图像,最大值为1:

其中,μx,σx和μy,σy分别代表重建图像和参考图像的像素值均值和标准差;Cov(X,Y)代表两图像间的协方差。

4.1 仿真实验分析

仿真实验使用经典Shepp-Logan人体头部模型来模拟和重建,并修改其模型参数为相应的相位参数来仿真纯相位头部模型［22］。

图2分别使用本文算法和FBP算法重建图像,其中,从左向右分别为180个、90个和36个投影角度下的重建图像;重建图像大小均为400×400像素;实验中迭代次数I,M分别为6和4。从图2可以看出,与FBP算法重建结果相比较,随着投影角度的减少,本文算法重建的图像质量下降不明显,依然可以保持较高的图像质量,且不产生显著的伪影。

图2 仿真实验的重建结果

将依次分析本文算法和FBP算法的定量指标并做相应对比。分别将相应算法在360个角度下的重建结果(满足数据完备性条件)作为相应算法的参考图像Y,结果如图3所示。

图3 仿真实验结果分析

由图3可以看出,随着投影角度的减少,本文算法在重建图像的信噪比、数值精确性以及与参考图像的相似度等方面比FBP算法有着明显的优势。在此基础之上,在稀疏投影角度情况下,本文算法可以较高精度地恢复出原始图像。

4.2 实际数据实验分析

实际实验投影数据来自日本高能加速器研究组织(High Energy Accelerator Research Organization)。同步辐射装置中X射线的能量是35 kV,检测器有450个。实验样品是一个有机玻璃体中包含13个圆柱状的通道,通道中分别填充水、甘油、乙醇。装置从上到下分别扫描500层,每层以0.5°的间隔扫描360次,共180°,每一层共采集360组数据。

本文中仿真实验使用本文算法和FBP算法分别

对最具代表性的第0层和第300层投影数据进行重建,采用等距间隔(分别为1°,2°,5°)获得180个、90个、36个角度进行重建。图4为实际数据实验重建图像,其中,从左至右依次为180个、90个和36个角度下的重建图像。从图4可以看出,在投影角度减少的情况下,本文算法在伪影控制方面比FBP算法拥有更好的效果,且可以得到不输FBP算法的图像重建效果。在衡量实际数据实验重建图像质量指标中,方便起见,本文只对第300层的结果进行分析并与FBP指标对比,结果如图5所示。

图4 实际数据实验重建图像

图5 实际数据实验结果分析

从图5可以看出,在实际数据实验中,较之FBP的分析结果,本文算法同样表现出更高的信噪比以及数值准确性,此外,与参考图像相似度方面也略胜一筹。证明本文算法在处理实际数据的图像重建问题中也有着不逊于FBP算法的性能。

5 结束语

基于压缩感知理论和恢复相位的策略,本文提出一种相衬CT迭代图像重建算法。该算法将折射角信息经预处理后,迭代重建出最终结果。在迭代重建部分,采用有序子集思想和自适应性权重,以增加对噪声的平滑作用和对不同投影数据的适应性。实验结果表明,对于欠采样的仿真数据和实际数据,本文算法有较好的图像重建质量和性能。由于目标函数带来的重建图像边缘锐度欠佳,因此后期需要改进目标函数构造来进一步优化图像质量。

［1]陈建文,高鸿奕,李儒新,等.X射线相衬成像［J］.物理学进展,2005,25(2):175-194.

［2]Davis T J,Gao D,Gureyev T E,et al.Phase-contrast Imaging of Weakly Absorbing Materials Using Hard X-rays［J］.Nature,1995,373(6515):595-598.

［3]Wilkins S W,Gureyev T E,Gao D,et al.Phase-contrast Imaging Using Polychromatic Hard X-rays［J］.Nature, 1996,384(6607):335-338.

［4]Momose A.DemonstrationofPhase-contrastX-ray Computed Tomography Using an X-ray Interferometer［J］.Nuclear Instruments and Methods in Physics Research SectionA:Accelerators,Spectrometers,Detectorsand Associated Equipment,1995,352(3):622-628.

［5]Chapman D,Thomlinson W,Johnston R E,etal.Diffraction Enhanced X-ray Imaging［J］.Physics in Medicine and Biology,1997,42(11).

［6]Snigirev A,SnigirevaI,KohnV,etal.Onthe Possibilities of X-ray Phase Contrast Microimaging by CoherentHigh-energySynchrotronRadiation［J］.Review ofScientificInstruments,1995,66(12): 5486-5492.

［7]Pfeiffer F,Bech M,Bunk O,et al.Hard-X-ray Dark-field Imaging Using a Grating Interferometer［J］.Nature Materials,2008,7(2):134-137.

［8]Wang Zhentian,Zhang Li,Huang Zhifeng,et al.An ART Iterative Reconstruction Algorithm for Computed Tomography of Diffraction Enhanced Imaging［J］.Chinese Physics C,2009,33(11):975-980.

［9]张 凯,朱佩平,黄万霞,等.代数迭代重建算法在折射衬度CT中的应用［J］.物理学报,2008,57(6): 3410-3418.

［10]李 镜,刘文杰,朱佩平,等.基于光栅相衬成像的扇束螺旋CT重建算法［J］.光学学报,2010,30(2): 421-427.

［11]Zhu Peiping,Zhang Kai,Wang Zhili,et al.Low-dose, Simple,and Fast Grating-based X-ray Phase-contrast Imaging［J］.Proceedings of the National Academy of Sciences,2010,107(31):13576-13581.

［12]李 镜,孙 怡.基于L1范数的微分相位衬度CT稀疏角度重建算法［J］.光学学报,2012,32(3):70-76.

［13]秦 峰,孙丰荣,宋尚玲,等.基于压缩感知的微分相衬CT迭代图像重建［J］.计算机应用,2013,33(6): 1732-1736.

［14]Candes E J,TerenceT.DecodingbyLinearProgramming［J］.IEEE Transactions on Information Theory, 2005,51(12):4203-4215.

［15]Candes E J,Romberg J K,Terence T.Stable Signal Recovery from Incomplete and Inaccurate Measurements［J］.Communications on Pure and Applied Mathematics,2006,59(8): 1207-1223.

［16]Candès E J,Romberg J,Terence T.Robust Uncertainty Principles:Exact Signal Reconstruction from Highly Incomplete Frequency Information［J］.IEEE Transactions on Information Theory,2006,52(2):489-509.

［17]Candès E J,Wakin M B.An Introduction to CompressiveSampling［J］.IEEESignalProcessing Magazine,2008,25(2):21-30.

［18]Sunaguchi N,Yuasa T,Huo Qingkai,et al.Convolution Reconstruction Algorithm for Refraction-contrast Computed Tomography Using a Laue-case Analyzer for Dark-field Imaging［J］.Optics letters,2011,36(3): 391-393.

［19]Hudson H M,Larkin R S.Accelerated Image Reconstruction Using Ordered Subsets of Projection Data［J］.IEEE Transactions on Medical Imaging,1994,13(4):601-609.

［20]De M B,Basu S.Distance-driven Projection and Backprojection［C］//Proceedings of IEEE Nuclear Science Symposium Conference.［S.l.］:IEEE Press,2002: 1477-1480.

［21]De M B,BasuS.Distance-drivenProjectionand Backprojection in Three Dimensions［J］.Physics in Medicine and Biology,2004,49(11):2463-2475.

［22]李涛涛,李 华,刁麓弘.复杂物体的X射线相衬成像三维仿真［J］.系统仿真学报,2009,21(1):32-35.

编辑 刘 冰

Sparse Projection Angle Phase Contrast CT Image Reconstruction Based on Refractive Angle

SI Kai1,SUN Fengrong1,SONG Shangling2,QIN Feng1
(1.School of Information Science and Engineering,Shandong University,Jinan 250100,China;
2.Department of Equipment,The Second Hospital of Shandong University,Jinan 250033,China)

It is not easy for conventional X-ray absorption contrast based Computed Tomography(CT)imaging technique to get a high quality image from the weak absorption materials comprising of light elements.X-ray phase contrast CT imaging is a CT imaging technique aiming at weak absorption materials and possessing ultrahigh resolution.However,X-ray phase contrast CT imaging takes too long time and requires too high X-ray radiation dose,which can set back the promotion of clinical applications.So it has significant meaning to study the image reconstruction for phase contrast CT using few projection views.This paper is under the theory of Compressed Sensing(CS),uses refractive angle information,aims to reduce the X-ray radiation dose,proposes a X-ray phase contrast CT image reconstruction algorithm.Experimental results show that compared with Filtered Back-projection(FBP)algorithm,this algorithm can use few projection views to obtain reconstructed images of higher quality,and can get higher Peak Signal to Noise Ratio(PSNR) and numerical accuracy in the actual data experiment.

X-ray phase contrast;Computed Tomography(CT);image reconstruction;Compressed Sensing(CS); refractive angle

司 凯,孙丰荣,宋尚玲,等.基于折射角的稀疏投影角度相衬CT图像重建［J］.计算机工程,2015, 41(3):262-268.

英文引用格式:Si Kai,Sun Fengrong,Song Shangling,et al.Sparse Projection Angle Phase Contrast CT Image Reconstruction Based on Refractive Angle［J］.Computer Engineering,2015,41(3):262-268.

1000-3428(2015)03-0262-07

:A

:TP391.41

10.3969/j.issn.1000-3428.2015.03.049

国家自然科学基金资助项目(61071053);山东省自然科学基金资助项目(ZR2014FM006)。

司 凯(1989-),男,硕士研究生,主研方向:医学图像处理;孙丰荣,教授;宋尚玲,博士;秦 峰,硕士研究生。

2014-03-03

:2014-05-18E-mail:sunfr@sdu.edu.cn