摘要:X射线三维成像技术是目前国内外X射线成像研究领域的一个研究热点. 但针对一些特殊成像目标, 传统X射线计算层析(CT)成像模式易出现投影信息缺失等问题, 影响CT重建的图像质量, 使得CT成像的应用受到一定的限制. 本文主要研究了基于光场成像理论的X射线三维立体成像技术. 首先从同步辐射光源模型出发, 对X射线光场成像进行建模; 然后, 基于光场成像数字重聚焦理论, 对成像目标场在深度方向上进行切片重建. 结果表明: 该方法可以实现对成像目标任一视角下任一深度的内部切片重建, 但是由于光学聚焦过程中的离焦现象, 会引入较为严重的背景噪声. 当对其原始数据进行滤波后, 再进行X射线光场重聚焦, 可以有效消除重建伪影, 提高图像的重建质量. 本研究既有算法理论意义, 又可应用于工业、医疗等较复杂目标的快速检测, 具有较大的应用价值.
关键词: X射线成像/
光场成像/
深度切片/
三维重建

English Abstract

全文HTML

--> --> -->

X射线发现至今, 各种Ｘ射线三维成像技术层出不穷, 尤其是同步辐射的广泛应用, 使得X射线层析成像在生物学、医药学、材料学等学科领域得到了广泛的应用^[1?4]. 传统Ｘ射线三维成像技术即X射线计算层析 (computer tomography, CT)成像技术^[5], 需要采集样品一定角度范围内的大量投影图, 然后经相应的CT重建算法重构, 再结合吸收衬度或相位衬度成像技术, 实现样品三维信息的精确重构^[6]. CT精确重建的基础是投影数据的完备性, 对于一些特殊场合或者样品, 无法采集完备的投影数据, 同时对样品进行大量的投影, 难以开展对时间分辨率要求较高的实时观测实验, 辐射剂量也相对较大, 使得CT成像技术的应用受到很大的限制. 另外, 若想获得特定视角下某一深度的所有图像信息, 必须首先通过完备的投影数据才能CT重建出所有的切片数据, 然后再通过其他三维重建软件, 方能依照深度信息对其进行进一步切片, 才能获得深度维度上的图像信息, 极大地增加了后期数据处理的复杂度和时效性.
为减少辐射剂量, 提高成像时间分辨率, 并准确获取样品的内部结构的深度信息, 上海光源^[7]、瑞士光源(PSI)^[8]和Spring-8光源^[9]通过不同的方式将同步辐射光折射后进而分为两束, 开展了基于双目视觉理论的X射线三维成像技术, 取得了一定的成果. 所谓双目视觉, 是基于同名像点视差原理, 由不同角度同时获取同一物体不同视角的数字图像, 根据投影几何三角法原理进行三维信息(深度信息)的恢复进而重构出物体三维几何信息^[10]. 但此方法成像之前需要对成像系统进行准确的标定, 在重构过程中, 需要对两个视角图像精确配准, 同时只有找到所有需要重构的目标在投影中的同名像素点, 才能进行重构. 对于一些复杂结构样品或者特征点不太明显的样品来说, 由此特征点匹配深度的重建算法准确度较差, 三维重建较为困难.
在可见光领域, 为了能快速准确地对目标场景进行三维成像, 光场成像技术得到迅速的发展^[11?13]. 在传统成像技术中, 成像系统采集到二维数字图像中每个像素都记录了到达该点的所有光线(或光子)的能量积分, 但并不记录这些光线的入射方向信息, 因此它只是三维场景(目标)的一个强度投影, 没有具体记录每一条光线的空间坐标和方向等, 丢失了三维结构的许多信息. 可见光领域的光场成像技术首先通过图像传感器的每个像素记录下来自任意不同角度的光线强度, 进而获得光场不同视角下的投影图像^[14], 即可通过传统的二维图像扩展记录包含有二维空间坐标强度和二维角度信息的四维光场信息; 然后通过计算成像重建技术及不同场景的重聚焦参数和相互遮挡效应, 对三维场景进行进一步的重建^[15?17].
本文以X射线传输特性为理论核心, 提出基于光场成像理论的X射线光场成像技术. 主要基于现有CT成像系统, 采集成像目标特定视角下的投影图像, 利用计算积分成像理论, 对空间中每一个体素进行重建, 并进一步在深度层面上对目标进行切片重建, 实现基于光场成像技术的样品三维结构的表征, 为X射线光场成像技术奠定理论基础并提供技术支持.

2

-->

2.1.光场重建原理

传统光场成像技术通过微透镜阵列^[12]或者相机阵列的方式^[18,19], 实现对目标场中每一条光线的二维强度及其二维入射方向信息进行采集, 即采集获得四维光场信息. 对于微透镜阵列式的光场相机, 一般依照图像传感器分辨率、微透镜阵列个数及二者之间的位置关系, 将成像镜头的入射光瞳划分为很多子孔径, 每一个子孔径坐标对应二维方向信息, 子孔径所成像代表此方向的二维投影信息. 在光场重建过程中, 将四维光场信息通过每一条光线与两个平行平面的交点坐标来进行参数化表征. 如图1所示, ${\varOmega _1}(u,v)$表示主透镜入射光瞳所在平面(即入射方向坐标平面), ${\varOmega _2}(s,t)$表示系统的像平面, 两个平面之间的距离为D, 光线$L(u,v,s,t)$的传播方向信息通过$(u,v)$表示, 空间坐标用$(s,t)$表示. 光场重建是将采集到的光场数据依靠两平面参数的对应关系重新投影到新的像平面上并进行积分. 在图1中, ${\varOmega _2}$面上的像等于${\varOmega _1}$与${\varOmega _2}$之间光场的积分, 即${\varOmega _2}$上每一点均为所有光线的经聚焦后积分结果. 新的像平面为${\varOmega _2}'$, 与${\varOmega _1}$面之间距离为D', 令$D' = \alpha D$, 对于每一条光线来说, 与${\varOmega _2}$和${\varOmega _2}'$都有交点, 由能量守恒定律可知:
图 1 光场成像原理示意图
Figure1. Schematic diagram of light field imaging principle.

$L'(u,v,s',t') = L(u,v,s,t).$

在新的像平面${\varOmega _2}'$上, 任一点$(s',t')$的像素值$I(s',t')$可以由所有经过该点的光线$L'(u,v,s',t')$对变量$u{\text{和}}v$积分获得, 并结合(1)式, 可得

$I(s',t') = \int {L'(u,v,s',t')} {\rm{ d}}u{\rm{d}}v = \int {L(u,v,s,t)} {\rm{ d}}u{\rm{d}}v.$

而又由几何对应关系$\displaystyle\frac{{s' - u}}{{D'}} = \displaystyle\frac{{s - u}}{D}$, 可得不同视角下, 经过${\varOmega _2}'(s')$的光线, 在${\varOmega _2}$上的坐标为

$s = \frac{{s'}}{\alpha } + u(1 - \frac{1}{\alpha }).$

同理可得

$t = \frac{{t'}}{\alpha } + v(1 - \frac{1}{\alpha }).$

将(3)式和(4)式代入(2)式可得新像面图像为

$\begin{split}I(s',t') = & \iint L\left( {u,v,u\left( {1 - \frac{1}{\alpha }} \right)} \right. + \frac{{s'}}{\alpha },\\& v\left. {\left( {1 - \frac{1}{\alpha }} \right) + \frac{{t'}}{\alpha }} \right){\rm{d}}u{\rm{d}}v.\end{split}$

对于成像系统来说, 只要确定了各个视角的投影信息及其所对应的角度坐标, 就可通过(5)式将图像投影到新的像面上, 实现对成像结果进行数字重聚焦.
2 -->

2.2.X射线光场成像模型

对于可见光成像, 来自成像目标的光线会辐射到空间的任意方向, 要想对成像目标进行成像, 必须先通过透镜对其进行聚焦, 进而通过微透镜阵列对其进行重采样^[20], 方可得到成像目标的多视角图像. 对于X射线光源来说, 发出的X射线传播方向一定, 在成像过程中, 也无需利用透镜对其进行聚焦, 只需让X射线穿透成像目标, 然后直接通过成像探测器采集图像, 以直接投影积分的方式实现对目标图像的获取^[21]. 亦即在一次成像结果中, 只有一个视角的成像信息, 不会出现多视角投影信息. 要想在一次成像结果中获得多视角投影信息, 可通过阵列光源在不同角度对成像目标进行照射, 然后通过对应的探测器采集相应视角的投影, 即达到多视角采集的目的. 虽然采用阵列光源形式一次成像可以获得更好的时间分辨率, 但是整个系统比较复杂且成本较高, 因此也可通过光源和探测器同步移动旋转或者只旋转成像目标的方式实现对成像目标的多视角投影图像采集. 图2为旋转成像目标实现多视角图像采集系统模型的截面示意图.
图 2 X射线光场成像系统模型示意图
Figure2. Schematic diagram of X-ray light field imaging system.

在图2所示的成像过程中, $(x,y,z)$为成像系统所在坐标系, $(x',y',z')$为成像目标所在坐标系, X射线自左向右沿着x轴传播, 两个坐标系原点均选在样品的旋转中心上. 由于X射线光场成像系统与传统光场相机在结构上完全不同, 所以四维光场信息的参数化表征方式也完全不同. 每一条光线可以表示为$L(x,y,z,\theta )$, x为深度信息, y和z为样品中的点在探测器像素平面上的投影坐标, 可以由CCD直接读取, $\theta $为投影视角信息. 对于基于同步辐射的平行X射线投影模型, 通过旋转样品的方式对样品进行多角度光场投影信息获取时, 光场信息角度采样只是在$xOy$平面内进行了一维采样. 对于X射线光场成像的“重聚焦”过程, 通过计算经过空间每一点的所有光线的累加积分实现. 在采集投影信息的初始时刻, 样品上某一点A坐标为$(x{'_1},y{'_1},z{'_1},\varphi )$, 经旋转$\theta $后, A移动到$A'(x{'_2},y{'_2},z{'_2},$$\varphi {\rm{ + }}\theta )$点的位置, 由几何对应关系, 经过$A{\rm{'}}$点的光线坐标$(x{'_2},y{'_2},z{'_2},\varphi {\rm{ + }}\theta )$与初始点A的光线坐标$(x{'_1},y{'_1},z{'_1},\varphi )$之间的关系为

$\left\{ \begin{aligned}&\left( \begin{aligned}&x{'_2}\\&y{'_2}\end{aligned} \right) = \left( \begin{aligned}&\;\;\;\,\cos \theta\;\;\;\sin \theta \\& - \sin \theta\;\;\;\cos \theta \end{aligned} \right)\left( \begin{aligned}&x{'_1}\\&y{'_1}\end{aligned} \right),\\&z{'_2}=z{'_1}.\end{aligned} \right.$

如图2所示, 当$\varphi = 0$时, 坐标系$(x',y',z')$和$(x,y,z)$重合. 根据(6)式, 在固定深度信息的情况下, 即$x'$为常数的情况下, 计算得出每一个投影角度下, 投影角度坐标$\theta $、投影图像坐标$(x,y,z)$以及成像目标中任一点空间$(x',y',z')$之间的关系为

$\left\{ \begin{aligned}&y = - x'\sin \theta {\rm{ + }}y'\cos \theta, \\&z = z'.\end{aligned} \right.$

通过光场积分理论可得, 点A$(x',y',z',0)$的像素值为

$I(x',y',z',0)=\int {L(x', - x'\sin \theta {\rm{ + }}y'\cos \theta ,z',\theta ){\rm{d}}\theta } .$

通过对每一点进行积分, 可以重建得样品空间每一个体素的值. 当固定$x'$时, 遍历$\theta $, 并按照(8)式计算积分将会得到如图2所示的$\varphi = 0$视角下, 固定深度$x'$处所对应的切片. 在重建过程中, 当初始坐标$(x',y',z')$与参考坐标系$(x,y,z)$不重合时, 即$x'$轴与$x$轴所在坐标系有一夹角$\varphi $, 首先需通过坐标变换, 将坐标系$(x',y',z')$变换到$(x,y,z)$坐标系, 然后积分遍历$\theta $将会得到$\varphi $度投影时的深度切片, 即此方法可以实现对样品不同视角的直接深度切片.

为了验证本文方法的准确性和正确性, 通过MATLAB模拟仿真的方式建立如图2所示的成像系统. 在系统中, 成像探测器的空间分辨率为256 × 50, 采用基于同步辐射的平行X射线束作为成像光源, 成像目标为堆叠了50层的MATLAB自带Shepp-Logan头骨模型, 在投影角度为0—180°范围内任意选取了64个角度, 各个角度及其投影数据如图3所示.
图 3 投影数据图　(a) 随机选取的64个投影角度; (b)图(a)中所有角度下的投影图中某一排像素图像组成的正弦图
Figure3. Projection data: (a) 64 random projection angles; (b) sinogram of some pixel image of projection in Fig.(a) under all angles.

图3(a)横坐标为投影编号, 纵坐标为投影所对应的角度坐标信息, 图3(b)为投影图像所对应的正弦图. 通过(8)式, 对投影数据进行重聚焦, 不同深度重聚焦结果和变视角的重聚焦结果如图4所示.
图 4 数字重聚焦结果　(a) Shepp-Logan模型; (a)中i线(b)、ii线(c)、iii线(d)和iv线(e)所在处的深度切片
Figure4. Digital refocus result: (a) Original Shepp-Logan phantom; depth slices where i line (b), ii line (c), iii line (d) and iv line (e) are located in Fig.(a).

光场成像重聚焦过程中, 通过选择不同的$x'$可以对样品进行深度切片. 如图4(b)为图4(a)中i线所在处的深度切片, 图4(c)为图4(a)中ii线所在处的深度切片, 其中在此深度切片中, 纵坐标为图像尺寸, 横坐标为成像目标堆叠的层数. 光场成像能够变换视角并在相应视角下对成像目标进行深度切片, 如图4(d)为图4(a)中iii线所在处的深度切片(视角$\varphi =45$°); 图4(e)为图4(a)中iv线所在处的深度切片(视角$\varphi = 90$°). 在光学成像的聚焦过程中, 超出景深范围的景物将会虚化, 使得准确聚焦景物更加突出, 成像效果更好, 当视场中的前后景物间存在相互遮挡时, 被遮挡的后景发出的光线将不会对成像有影响, 使得其对前景的重聚焦图像没有任何贡献. 但是X射线成像属于完全透射的成像模式, 如图2所示, 经过前景的光线必然会携带后景的信息. 在重聚焦过程中, 所在重聚焦深度的物体会被聚焦, 其他背景信息会处于离焦状态, 过多的背景信息会淹没图像中真正的有用信息, 使得切片中整个图像, 尤其是图像的边缘不再清晰, 对成像非常不利, 如图4(b)—(e)的重聚焦结果中, 虽然能从深度切片中看到样品的特征信息, 但是Shepp-Logan模型特征信息非常不清晰, 边界也很模糊. 尤其对于图4(e)的切片, 原始数据的细节信息很丰富, 但是此切片中很模糊, 细节信息未得到很好的体现. 为提高成像效果, 增加重聚焦图像及其边缘清晰程度, 通过对原始数据进行滤波的方式, 可以消除背景的影响, 图5为添加了R-L滤波器后光场重聚焦结果.
图 5 经R-L滤波器滤波后的数字重聚焦结果　(a) Shepp-Logan模型; (a)中i线(b)、ii线(c)、iii线(d)和iv线(e)所在处的深度切片
Figure5. Digital refocus result after filtering by R-L filter: (a) Original Shepp-Logan phantom; depth slices where i line (b), ii line (c), iii line (d) and iv line (e) are located in Fig.(a).

经滤波后的原始数据再进行数字重聚焦, 可以得到任意视角和任意深度下的样品内部切片, 如图5(b)和图5(c)所示的0°视角情况下不同深度的样品内部信息, 也可以得到如图5(d)和图5(e)所示的$\varphi = 45$°和$\varphi = 90$°等不同视角的深度图像. 在各个切片图像中, 很好地恢复出Shepp-Logan模型图像的特征, 图像的清晰度和边缘得到了极大提升, 如图5(e)可以很好地展示出iv线所在处样品的丰富信息. 另外, 也可以先对样品进行重聚焦重建, 根据重聚焦参数$x'$确定重聚焦结果中特征信息在样品中的深度. 在取得所有图像的切片之后, 同样可以通过相关三维重建软件, 对切片进行三维立体结构的重建进一步恢复得到样品的整体三维结构.

本文利用同步辐射光源的准平行光模型, 基于光场成像理论系统地研究了完全透射成像模型中成像目标的光场重建算法; 在获取任意角度目标场投影信息的前提下, 通过几何追迹的方式获取经过成像目标中任一体素的光线, 利用积分成像理论, 恢复得到该体素的图像灰度值, 进而得到任一深度位置上的切片图像, 并进一步研究了光学成像中离焦景物所引入背景噪声的去除等问题.
理论分析和仿真实验结果表明, 本文在获得样品的深度切片信息的基础上, 同时也可以对典型目标进行深度测量. 将原始数据滤波以后, 获得不受背景影响的高信噪比样品的内部结构. 相对于传统解析层析重建算法的等角度扫描重建, 本文所示的方法可以任意扫描方式获取原始数据, 并可直接获取成像目标任一深度切片, 对于一些复杂成像目标可以只对感兴趣区域进行扫描进而恢复得到深度信息, 具有更强的灵活性和适应性. 本文通过旋转样品的方式采集光场数据, 但是此方式成像时间分辨率较低, 对于快速成像的场合, 利用光场成像数字重聚焦理论, 可以使用阵列X射线光源照射同时采用阵列探测器采集成像数据, 经一次同步曝光, 校准后的系统同样可实现目标场的光场三维重建, 具有很高的时间分辨率.