基于迭代重建算法的X射线光栅相位CT成像∗

戚俊成 陈荣昌 刘宾 陈平 杜国浩 肖体乔

1)(中北大学信息与通信工程学院,太原 030051)

2)(中国科学院上海应用物理研究所,上海 201800)

基于迭代重建算法的X射线光栅相位CT成像∗

戚俊成1)†陈荣昌2)刘宾1)陈平1)杜国浩2)肖体乔2)‡

1)(中北大学信息与通信工程学院,太原 030051)

2)(中国科学院上海应用物理研究所,上海 201800)

(2016年8月18日收到;2016年12月5日收到修改稿)

基于光栅干涉仪的X射线成像技术可以同时获得样品内部的吸收信息、相位信息和散射信息,既保持了传统X射线衰减成像的优点,又拥有相衬成像和散射成像的优势.然而基于传统CT重建算法的X射线光栅成像需要采集大量完整的原始投影数据,数据采集时间过长从而使得物体接受很大的辐射剂量,难以在实际中应用.提出基于传统代数迭代重建算法的光栅成像技术.该方法利用现有X射线光栅成像系统采集少量原始投影数据,基于传统代数迭代重建算法,对旋转变化的相位数据进行CT重构,同时基于傅里叶变换的方法对微分相位数据进行相位恢复.模拟和实验结果表明,基于少量或不完备的原始投影数据,该方法能够准确重构成像对象的吸收、相位和散射三维信息,同时还能对微分相位切片进行高信噪比的相位恢复,得到样品折射率实部衰减率,为X射线光栅成像技术在工业、生物和医学诊断等领域的应用提供理论和技术支撑.

X射线光栅成像,迭代重建算法,相位恢复

1 引 言

自从X射线发现至今,X射线成像一直受到广泛的关注.X射线光栅成像技术对X射线源的相干性要求较低,并可以同时获得样品的吸收、散射和相位信息[1].更重要的是,基于光栅的X射线相衬成像方法不受视场限制,可以实现大视场成像.因此,该方法被视为最有潜力从实验室走向实际应用的X射线相衬成像方法,受到了广泛的的关注,也具有广阔的实际应用前景.

传统光栅成像技术的重建算法是以希尔伯特变换为核心的反投影算法,是一种解析重建算法,要求具有完整的投影数据.在数据采集过程中,以步进扫描的方式采集数据需要花费大量时间采集大量投影作为原始数据,使得样品接收大量的辐射.针对这个问题,国内外科学家采取了很多办法来减少投影数,从而降低样品所受辐射剂量.Momose等[2]和Wen等[3]采用一次曝光的方法采集数据,即样品每旋转一个角度,采集一张叠栅条纹图像,然后将叠栅条纹图像进行傅里叶变换,对比有样品和没有样品情况下经傅里叶变换后的各个参数,就可以得到样品的各种信息.这种方法的优点是不需要相位步进扫描,时间分辨率较高,但是成像的空间分辨率受到叠栅条纹周期限制,并且对光栅和探测器的要求也很高,要求相位光栅和吸收光栅完全匹配.Zhu等[4]采用反投影的方法,在样品0◦—360◦旋转的同时采集数据,然后对比样品旋转角度差为180◦时的两幅图像,从而得到样品信息.这种方法在CT过程中省去了步进扫描,使得实验效率大大提升,但是无法获得散射信息.Jensen等[5]采用交错光栅步进的方法采集数据,在样品旋转的同时进行光栅扫描.Zanette等[6,7]提出了交错相位步进扫描法,即在样品旋转的同时进行相位步进扫描,从多个样品旋转角度中获取样品的信息.在同等辐射剂量下,此方法可以获得较好的结果,但是也存在投影角度不精确的问题,当投影数较少时,远离旋转中心的样品高频信息会丢失,而且此方法较难实现.

此外,光栅成像在实际应用中经常无法获得完整的投影数据,主要原因是辐射剂量、经济成本、无法检测或很难检测到完全的投影数据等.对于很大的样品,要求探测器和光栅面积足够大才能重建,难以实现;有时投影不是均匀地分布在180◦或360◦范围内,这时滤波反投影(filtered back projection,FBP)方法受到了限制.

综上所述,本文针对X射线光栅成像技术存在的问题,通过传统代数迭代重建算法(algebraic reconstruction technique,ART),利用少量或者非完整的投影数据,对快速、低剂量光栅成像技术的CT重建和相位恢复算法展开了研究.利用传统ART对由光栅成像所获得的吸收、散射和旋转变化的一阶微分相位信息进行了高精度重建.对一阶微分相位切片进行了高信噪比的相位恢复,为X射线光栅成像技术在生物医学、材料和工业等领域的研究应用提供技术支持.

2 原 理

2.1 光栅成像技术

基于准平行光束的X射线光栅成像系统主要由相位光栅G1、吸收光栅G2和探测器组成(图1)[8].当X射线光束穿过相位光栅G1后,由于泰伯效应(光栅自成像效应)[9],在G1后的泰伯距离处,将会得到光栅G1的自成像.在此像面上放置一个周期与自成像条纹对应的吸收光栅G2,即可得到叠栅干涉条纹.当X射线透过放置于相位光栅前的样品后,X射线被物体吸收、折射和散射,相位光栅的入射X射线波前发生变化,从而导致相位光栅的泰伯像条纹发生强度变化和形状扭曲.光栅成像系统的作用之一就是定量测出此泰伯像条纹的扭曲度和光强度的变化,最终得到X射线穿过样品的吸收、相位和散射信息.

为了从探测器所得到的信息中分离出相位信息和其他信息,采用相位步进技术使一块光栅以步进的方式扫描另外一块光栅,获得每一个像素点上随着步进数变化的光强位移曲线[10].在忽略高次项的情况下,位移曲线方程可以傅里叶变换展开为[1,11,12]

式中(m,n)为像素点的坐标,k=2π/d,d为光栅G2的周期,yg为相位步进扫描过程中相位光栅所处的位置,ai为振幅系数,ϕi为相对应的相位系数.

图1 X射线光栅成像系统基本布局图Fig.1.Layout for the X-ray grating imaging system.

样品的吸收信息与每一像素点位移曲线的0阶振幅系数a0有关;微分相位信息由位移曲线的相位系数ϕ1表征[11];散射信息主要反映了样品的局部散射能力,样品的散射信息可由条纹的可见度V=a1/a0表征[1].

2.2 光栅成像的物理模型分析

X射线光栅成像技术的一大优点就是通过一组原始数据可以同时获得样品的吸收、相位和散射信息,而这三种信息的表征形式是不同的.吸收信息和散射信息的表征形式是样品的吸收系数或者散射因子的线积分,而在整个CT扫描过程中,吸收系数的整体分布随着样品一起转动,即吸收系数都是旋转不变的;对于大部分成像样品来讲,也可以认为散射因子是旋转不变的或者近似旋转不变的.所以,对于吸收信息和旋转不变的散射信息,只需要通过实验获得相应的CT投影数据,利用传统迭代重建算法就可以对其进行CT重建.

相位信息的数据采集过程如图2所示,o′x′y′z′为样品所在坐标系,随着样品一起转动,oxyz是光栅和X射线光源所在坐标系,是静止不动的.在CT数据采集过程中,样品折射率实部衰减率δ的整体分布随着样品一起转动,即δ的分布是旋转不变的.但是光栅所在坐标系相对于样品坐标系不是静止的,所以δ沿光栅坐标系y方向的一阶偏导的分布是旋转变化的.因此,微分相位信息不能像吸收信息或者散射信息一样,直接应用传统迭代CT重建算法进行CT重构.

图2 相位信息数据采集过程示意图Fig.2.Schematic of phase acquisition.

通过对X射线在样品中传播过程的分析可知,样品的相位信息取决于X射线穿过样品后的折射角,而折射角为∂δ/∂y沿着X射线传播方向的积分.在CT数据采集过程中,随着投影角度ϕ的变化,投影数据、投影角度和δ一阶偏导数之间的关系为[13]

代表投影角度为ϕ时恢复得到的相位导数和分别代表初始δ分布沿着垂直于X射线的传播方向和平行于X射线的传播方向上的一阶偏导数.通过(2)式可以将旋转变化的的投影信息转变为旋转不变的的投影信息,即将旋转变化的微分相位信息转变成旋转不变的微分相位信息.和

2.3 代数迭代CT重建

对于滤波反投影CT重建算法而言,需要完整并且角度分布均匀的投影数据,而迭代重建适合投影数据不完全或投影间隔不均匀等场合的图像重建[14].传统代数迭代CT重建方法的基本思想是:首先将图像离散地分解为N=n×m的矩形区域,每一个区域代表一个像素,同时给定一个初始图像估计值f(0),然后计算f(0)的投影数据,并用该投影数据与实际投影数据的差来校正图像,得到一次近似图像f(1).按照同样的方法再由f(1)求得二次近似图像f(2),当投影到最后,得到f(N),即完成一次迭代过程.第二次迭代时,以f(N)为初始值并再次投影,如此继续,直到满足预定条件为止,其迭代公式为

式中b为迭代次数,aij为投影矩阵系数,pi为投影数据向量,xj为重建后的图像的第j个像素值,λ为松弛因子,取值范围为(0,2).

2.4 相位恢复

利用迭代算法重建得到的微分切片数据,构建新的二维函数f(y′,z′)[15]:

由傅里叶变换的微分特性可得

即[16]

式 中F[f(y′,z′)](k,l)是f(y′,z′)的 傅 里 叶 变 换,(k,l)是傅里叶变换后的坐标,因此相位切片数据可以表示为

式中F和F−1分别为傅里叶变换和傅里叶逆变换.因此,利用δ沿着y′和z′方向的一阶偏导数切片数据可以准确恢复得到δ.

综合以上分析,通过ART可实现光栅成像技术与迭代重建算法的完美结合,并对δ进行高信噪比的相位恢复.

3 结果分析

3.1 模拟结果

在对光栅成像数据进行CT重构时,最需要解决的问题是微分相位数据的CT重建和相位信息的恢复.为了验证本文算法的准确性,对理想相位模型进行光栅成像模拟并对其进行相应重构和恢复.改进Shepp-Logan模型,得到理想相位模型(图3(a)和图3(b)).基于同步辐射X射线光源模型,按照表1所示的参数,模拟光栅成像过程及其相应的相位步进扫描,得到每一个像素点有样品和无样品时的位移曲线,通过光栅成像信息恢复算法公式((1)式),得到折射率实部衰减率导数∂δ/∂y的投影正弦图(图3(c)).通过(2)式得到投影正弦图在垂直和平行于X射线传播方向上的微分投影信息,即和的投影信息,利用(3)式对其进行CT重构后,得到和的切片(图4). 最后结合两个方向上的一阶微分信息进行相位恢复,得到图5所示相位恢复结果.

表1 模拟参数Table 1.Simulation parameters.

图3 (a)Shepp-Logan模型;(b)(a)中白线位置处的δ分布图;(c)微分相位正弦图Fig.3.(a)Shepp-Logan phantom;(b)the profile at the position marked with white line in(a);(c)the sinogram of the differential phase information.

图4 微分相位信息重建结果 (a)δ沿着y′方向微分切片;(b)δ沿着z′方向微分切片Fig.4.Reconstruction result with iterative reconstruction algorithm:(a)Differential phase image of phantom along y′direction;(b)differential phase image of phantom along z′direction.

图5 折射率迭代算法重建结果 (a)δ切片图;(b)(a)中白线位置处的δ分布图Fig.5.Reconstruction result of the refractive index with iterative reconstruction algorithm:(a)The slice of δ;(b)the profile at the position marked with white line in(a).

图6 折射率滤波反投影重建结果 (a)δ切片图;(b)(a)中白线位置处的δ分布图Fig.6.Reconstruction result of the refractive index with FBP:(a)The slice of δ;(b)the profile at the position marked with white line in(a).

利用旋转不变的折射率实部衰减率的一阶偏导数投影信息,通过传统迭代重建算法,可以得到图4所示的微分相位信息的切片数据.由图4(a)和图4(b)所示切片信息可以很明显地看到样品的边界信息,且求导方向互相垂直.最后恢复得到的δ信息切片如图5所示.滤波反投影的重建结果如图6所示.对比图5(a)和图6(a)所示切片可以发现,由于投影数量较少,两种方法均会出现一定量的伪影,均可依据图像灰度区分不同δ所对应的区域,但迭代重建算法的图像质量明显高于滤波反投影的重建结果.从图5(b)和图6(b)所示对应切片中相同位置处的δ分布图可以看出,迭代重建算法得到切片的定量信息更加准确,信噪比更高,具有重要的实际应用意义.

3.2 实验结果

实验在上海光源(SSRF)X射线成像及生物医学应用光束线站(BL13W1)上进行,该线站采用1.9T的Wiggler光源和液氮冷却双晶单色器,能提供8—70 keV的单色X射线[17].成像系统由旋转样品台、相位光栅、吸收光栅、纳米位移台和探测器组成,如图7所示.光栅成像系统距离光源点约34 m,在此处X射线光束基本可看作准平行光,具有较好的空间相干性[18].相位光栅安装在PI公司的P-611.ZS型纳米精度压电陶瓷位移平台上,用以实现纳米精度相位步进扫描.吸收光栅安装在一个精密转台上,可根据实验情况实时同步调整两光栅刻线之间的夹角以改变叠栅条纹的周期.两光栅之间的水平距离和垂直高度可分别通过安装在其底座上的精密位移台远程自动控制.在采集图像过程中,使用Photonic Science公司生产的9µm X射线CCD记录叠栅条纹变化.

为突出验证光栅成像和迭代重建算法的优势,实验中采用能量为20 keV的X射线对样品在0◦—180◦范围内作CT旋转扫描,CT扫描间隔为3◦;移动吸收光栅,在每一个CT投影角度利用多步相移法对相位光栅进行相位步进扫描,相位步进扫描的步长为0.6µm,即每周期扫描4次,具体实验参数如表2所示,成像样品为尼龙6(PA6)包裹的聚甲基丙烯酸甲酯(PMMA)棒材,样品在20 keV光束能量时,折射率等信息如表3所示.

图7 上海光源光栅成像装置图Fig.7.Grating based X–ray imaging system at SSRF.

通过多步相移法,利用光栅成像系统采集得到原始投影信息,按照(1)式所示的方法进行信息提取后,可以分别得到样品的吸收、相位和散射信息.吸收和散射信息直接通过传统代数迭代重构算法重构后,可得到样品的吸收和散射信息切片(图8);将旋转变化的微分相位信息进行变换后,利用迭代重建算法重建,可得到δ的微分切片如图9(a)和图9(b)所示.再对图9(a)和图9(b)所示微分切片进行相位恢复,最后得到图9(c)所示的样品折射率实部衰减信息.

表3 样品信息列表Table 3.Information of samples.

图8 重建结果 (a)吸收切片;(b)散射切片Fig.8.Experimental results of reconstruction:(a)Slice of absorption information;(b)slice of scattering information.

图9 相位信息重建结果 (a)δ沿着y′方向微分切片;(b)δ沿着z′方向微分切片;(c)恢复后的δ信息Fig.9.Reconstructed results of phase information:(a)Differential phase image of the sample along y′direction;(b)differential phase image of the sample along z′direction;(c)the slice after phase retrieval.

由以上重构结果可以看出,在少量投影的情况下,通过传统代数迭代重建算法可以获得样品各种衬度的准确信息.在图8(a)所示的吸收切片中,样品为低Z材料,吸收系数很小,只能看到由于同步辐射X射线成像边缘效应而形成的边界信息,样品的内部信息与空气很难通过图像予以区别.图8(b)所示的散射信息中只有各种样品的边界信息,这是因为样品是均匀材料,只有在不同材料的交界处才会对X射线形成散射.因此,通过吸收和散射信息很难将PA6,PMMA和空气等区分开来.而对于相位信息,虽然两种材料的δ差别很小(表3),但是相位恢复后由图9(c)可准确区分PA6和PMMA两种高分子材料.实验样品的成分和结构更加复杂时,不同衬度的切片将会突出展示样品的不同信息,使得成像结果更加丰富,最终形成吸收、相位、散射三位一体的成像模式.

4 结 论

本文基于上海光源同步辐射平行光模型,利用传统的代数迭代重建算法系统地研究了少量投影情况下X射线光栅成像的信息重建算法,重点研究了一阶微分相位信息的CT重建和样品折射率实部衰减率δ的恢复问题.通过变换的方式,将随着投影角度变化而变化的一阶微分相位信息转变为旋转不变的信息,即沿着光线和垂直于光线两个方向的微分信息,通过传统的ART,对其直接进行重建得到了δ的一阶微分切片,最后再通过傅里叶变换恢复得到δ的实际分布信息.

理论分析和实验结果表明,该方法不仅可以获得样品的吸收切片、沿着两个不同方向的微分相位切片和散射切片数据,还可以获得高信噪比的相位数据,在X射线光栅成像中具有可行性.相对于传统的反投影相位重建算法,基于ART的相位重建算法在少量投影情况下,在光栅成像物理模型和降低辐射剂量等方面都有所改善;与传统积分对一阶微分信息的恢复算法相比,该方法不需要考虑微分数据的边界条件,也没有误差累积效应,可以恢复得到信噪比非常好的相位数据;由于相位恢复时,使用的两个方向的微分数据均为同一次实验结果,也不存在图像配准等问题.因此,本文提出的快速光栅成像方法既保留了光栅成像和ART重建算法的优点,又可以实现在缺角度或者光束不完全覆盖样品等条件下的X射线光栅CT成像,在生物医学、材料和工业等要求低辐射剂量、快速和高精度的领域具有很好的应用前景.

感谢上海光源BL13W1线站工作人员在实验工作中的帮助.

[1]Pfeiffer F,Bech M,Bunk O,Kraft P,Eikenberry E F,Brönnimann C,Grünzweig C,David C 2008Nat.Mat.7 134

[2]Momose A,Yashiro W,Maikusa H,Takeda Y 2009Opt.Express17 12540

[3]Wen H H,Bennett E E,Kopace R,Stein A F,Pai V 2010Opt.Express35 1932

[4]Zhu P P,Zhang K,Wang Z L,Liu Y J,Liu X S,Wu Z Y,McDonald S A,Marone F,Stampanoni M 2010Proc.Natl.Acad.Sci.107 13576

[5]Jensen T H,Bech M,Zanette I,Weitkamp T,David C,Deyhle H,Rutishauser S,Reznikova E,Mohr J,Feidenhans’l R,Pfeiffer F 2010Phys.Rev.B82 214103

[6]Zanette I,Bech M,Pfeiffer F,Weitkamp T 2011Appl.Phys.Lett.98 094101

[7]Zanette I,Bech M,Rack A,Le Duc G,Tafforeau P,David C,Mohr J,Pfeiffer F,Weitkamp T 2012Proc.Natl.Acad.Sci.109 10199

[8]Chen B,Zhu P P,Liu Y J,Wang J Y,Yuan Q X,Huang W X,Ming H,Wu Z Y 2008Acta Phys.Sin.57 1576(in Chinese)[陈博,朱佩平,刘宜晋,王寯越,袁清习,黄万霞,明海,吴自玉2008物理学报57 1576]

[9]Talbot H F 1936Philos.Mag9 401

[10]Qi J C,Ren Y Q,Du G H,Chen R C,Wang Y D,He Y,Xiao T Q 2013Acta Opt.Sin.33 1034001(in Chinese)[戚俊成,任玉琦,杜国浩,陈荣昌,王玉丹,和友,肖体乔2013光学学报33 1034001]

[11]Bech M,Jensen T H,Bunk O,Donath T,David C,Weitkamp T,Le Duc G,Bravin A,Cloetens P,Pfeiffer F 2010Zeitschrift Fur Medizinische Physik20 7

[12]Momose A,Kawamoto S,Koyama I,Suzuki Y 2004Developments in X-Ray Tomography IV5535 352

[13]Zhu P P,Wang J Y,Yuan Q X,Huang W X,Shu H,Gao B,Hu T D,Wu Z Y 2005Appl.Phys.Lett.87 264101

[14]Yang F Q,Zhang D H,Huang K D Wang K,Xu Z 2014Acta Phys.Sin.63 058701(in Chinese)[杨富强,张定华,黄魁东,王鹍,徐哲2014物理学报63 058701]

[15]Kottler C,David C,Pfeiffer F,Bunk O 2007Opt.Express15 1175

[16]Arnison M R,Larkin K G,Sheppard C J R,Smith N I,Cogswell C J 2004J Microsc.214 7

[17]Xiao T Q,Xie H L,Deng B,Du G H,Chen R C 2014Acta Opt.Sin.34 0100001(in Chinese)[肖体乔,谢红兰,邓彪,杜国浩,陈荣昌2014光学学报34 0100001]

[18]Qi J C,Ye L L,Chen R C,Xie H L,Ren Y Q,Du G H,Deng B,Xiao T Q 2014Acta Phys.Sin.63 104202(in Chinese)[戚俊成,叶琳琳,陈荣昌,谢红兰,任玉琦,杜国浩,邓彪,肖体乔2014物理学报63 104202]

PACS:42.25.Kb,41.60.Ap,07.60.Ly,87.59.–e DOI:10.7498/aps.66.054202

Grating based X-ray phase contrast CT imaging with iterative reconstruction algorithm∗

Qi Jun-Cheng1)†Chen Rong-Chang2)Liu Bin1)Chen Ping1)Du Guo-Hao2)Xiao Ti-Qiao2)‡

1)(School of Information and Communication Engineering,North University of China,Taiyuan 030051,China)
2)(Shanghai Institute of Applied Physics,Chinese Academy of Sciences,Shanghai 201800,China)

18 August 2016;revised manuscript

5 December 2016)

Grating based X-ray imaging technology is a coherent imaging technique that bears tremendous potential in threedimensional tomographic imaging of weak absorption contrast specimens.Three kinds of contrast information including absorption,phase and scattering can be retrieved separately based on a single set of raw projections.However,the grating based X-ray imaging with the conventional phase-retrieval method using the conventional phase-stepping approach andfiltered back projection(FBP)reconstruction algorithm require large amounts of raw data,so that long exposure time and large amounts of radiation dose is accepted by the sample.According to the traditional grating based X-ray imaging system,we propose a low dose,fast,multi-contrast CT reconstruction approach based on the iterative reconstruction algorithm that optimizes dose efficiency but does not share the main limitations of other reported methods.Prior to reconstruction,firstly,the projections are acquired with the phase stepping approach and multi-contrast projections are retrieved from the raw data by conventional retrieval algorithm.Then the rotational variable differential phase projections are converted to rotational invariable projections by means of decomposing the differential phase projections into the rotational invariable projections in two mutually perpendicular derivative directions via the transformation of coordinates.Finally,the absorption,phase and scattering information are reconstructed with the iterative reconstruction algorithm and the phase is retrieved based on the fast Fourier transform(FFT).We validated and assessed the phase reconstruction approach with a numerical simulation on a phase Shepp-Logan phantom.The experiment was performed at the X-ray imaging and biomedical application beam line(BL-13W)in the Shanghai Synchrotron Radiation Facility(SSRF)where 20 keV X-ray from a Si(111)monochromator is emitted.The X-ray interferometer was positioned at 34 m from the Wiggler source.The images were recorded with a scintillator/lens-coupled CCD camera with 2048 pixel×2048 pixel resolution and an effective pixel size of 9µm.The numerical tests and the experimental results demonstrate that,for the small radiation dose deposited in the sample,the iterative reconstruction algorithm provides phase reconstructions of better quality and higher signal to noise ratio than the conventional FBP reconstruction algorithm,and also provides the multi-contrast 3D images,including absorption image,phase image and scattering image.This development is of particular interest for applications where the samples need inspecting under low dose and high speed conditions,and will play an important role in the nondestructive and quantitative imaging in the industry,biomedical and medical diagnosisfields.

grating based X-ray imaging,iterative reconstruction algorithm,phase retrieval

PACS:42.25.Kb,41.60.Ap,07.60.Ly,87.59.–e

10.7498/aps.66.054202

∗国家自然科学基金(批准号:11375257,61301259,U1232205)、中北大学校学科研究基金(批准号:2015110246)和山西省自然科学基金(批准号:2015021099)资助的课题.

†通信作者.E-mail:qijuncheng@nuc.edu.cn

‡通信作者.E-mail:tqxiao@sinap.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.11375257,61301259,U1232205),the Foundation of North University of China(Grant No.2015110246),and the Natural Science Foundation of Shanxi Province,China(Grant No.2015021099).

†Corresponding author.E-mail:qijuncheng@nuc.edu.cn

‡Corresponding author.E-mail:tqxiao@sinap.ac.cn