CN101467889B - 光栅剪切相位衬度ct成像数据采集和重建方法 - Google Patents

Abstract

一种光栅剪切相位衬度CT成像数据采集和重建方法，利用分析光栅相对于相位光栅的剪切技术，测得光强随分析光栅剪切位移变化的位移曲线，在分析光栅处于位移曲线中值点位置(腰位)、样品转轴与光栅栅条平行时，采集一套样品旋转360度的CT投影数据；从360度的CT投影数据中选取方向相反的投影光线数据进行配对组合，分离出样品旋转180度的吸收投影数据和折射投影数据，再重排变换，形成样品旋转180度的平行束吸收投影数据和折射投影数据；根据旋转不变性要求对折射投影数据进行加权，得多套满足旋转不变性要求的180度平行束投影数据，再依雷当逆变换重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

Description

光栅剪切相位衬度CT成像数据采集和重建方法

技术领域

本发明涉及计算机体层成像技术领域，特别是一种光栅剪切相位衬度CT成像数据采集和重建方法。

背景技术

传统的计算机体层成像(Computed Tomography，简称CT)是现代科学技术对人类健康的巨大贡献。经过三十多年的发展，该项技术在数据采集、重建方法和应用等多方面都取得了令人瞩目的进展，成为医学、生物、考古、材料科学、安全检查和工业无损检测等领域不可或缺的重要工具。传统CT原理的物理基础是物质对X射线吸收的差别，数学基础是雷当(Radon)变换。因此，传统CT技术是X射线吸收衬度旋转投影成像技术与雷当变换的结合。

基于吸收机制的X射线成像技术仅对人体骨骼观察的比较清楚，而对人体软组织成像模糊，特别是难以观察到人类体内毫米量级的早期恶性肿瘤(癌症)。这其中的原因在于人体软组织主要是由轻元素组成的，轻元素对于硬X射线几乎没有吸收，就象可见光透过水中的一个小玻璃球，几乎没有留下可以察觉的痕迹。目前癌症已经成为人类的主要杀手，实现癌症的早期诊断，是人类当前战胜这一恶魔的关键。因而需要发展一种方法，能清楚地分辨癌组织和正常组织，这就好象要想办法看清楚水中的一个小玻璃球。科学家发现，轻元素物质引起X射线相位改变的幅度是其对X射线吸收值的一千倍到十万倍，利用相位信号发展的X射线相位衬度成像，特别适合观察轻元素构成的物体。

近年来发展起来的X射线相位衬度投影成像方法，利用相位衬度比吸收衬度高得多的特殊性质，对生物软组织和轻元素样品获得了衬度非常高的成像结果。由于X射线相位衬度投影成像方法具有的巨大优势和广阔的发展前景，所以在国际上得到了迅速发展。目前已经发展了四种相位衬度投影成像方法，它们分别是晶体干涉仪成像方法、相位传播成像方法、衍射增强成像方法和光栅剪切成像方法。在四种相位衬度投影成像方法中，虽然光栅剪切成像方法是本世纪才发展起来的，但是它却是最有发展前途的投影成像方法。它不但能利用角度分辨能力获得与折射角成正比的图像衬度，而且具有和常规X射线光源结合、向医学临床应用发展的前景。因此，将光栅剪切X射线相位衬度旋转投影成像技术与雷当变换结合就成为人们进行研究的课题。

目前国际上已经提出的光栅剪切相位衬度CT成像方法，利用X射线穿过样品产生的折射角(也是相位梯度)信号，已经对生物软组织成像获得了非常高的衬度。然而，目前的光栅剪切相位衬度CT投影数据的采集方法和重建方法太麻烦。例如，需要将分析光栅分别固定在位移曲线的五个不同位置上，依次采集五套样品旋转180度的投影数据，也就是说在数据采集过程中，分析光栅需要移动五次，样品至少要依次旋转5个180度，才可能从采集的投影数据中解出样品的相位信息。由此可见，这种投影数据的采集方法比传统CT方法复杂。除此之外，已经提出的重建方法只能重建出样品的吸收系数和折射率。显然，目前的光栅剪切相位衬度CT投影数据的采集方法和重建方法不适合推广到人体成像。

发明内容

本发明的目的是克服现有技术存在的缺陷，公开一种光栅剪切相位衬度CT成像数据采集和重建方法，操控简便，投影数据采集快，不仅可以重建出样品的吸收系数和折射率，而且可以重建出衬度更高的折射率一阶导数和折射率二阶导数。因此，本发明具有将相位衬度三维成像技术推广到医学临床的前景。

为达到上述目的，本发明的技术解决方案是：

一种光栅剪切相位衬度CT成像数据采集和重建方法，其在分析光栅处于位移曲线中值点位置(腰位)，样品转轴与光栅栅条平行的条件下，采集一套样品旋转360度的CT投影数据；从360度的CT投影数据中选取方向相反的投影光线数据进行配对组合，分离出样品旋转180度的吸收投影数据和折射投影数据，再经重排变换，形成样品旋转180度的平行束吸收投影数据和折射投影数据；根据旋转不变性要求对折射投影数据进行加权，获得多套满足旋转不变性要求的180度平行束投影数据，再根据雷当逆变换重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

所述的光栅剪切相位衬度CT成像数据采集和重建方法，其具体步骤如下：

步骤S1，准备光栅剪切成像条件：

用相干的X射线光束照射相位光栅，根据光学中的泰保效应，在相位光栅后面会产生自成像，自成像的空间频率是相位光栅的空间频率的一倍或者二倍，在自成像位置，放置空间频率和自成像空间频率相同的分析光栅，沿着垂直于光束和光栅栅条的方向，移动分析光栅，同时利用光强探测器测得光强随着分析光栅位置移动而起伏的位移曲线，并拍摄无样品的背景光分布；

步骤S2，采集样品旋转360度的相位衬度CT投影数据：

将分析光栅位置调整在位移曲线的中值点(腰位)，在相位光栅前安装样品转台，调整样品转台旋转轴方向，使其平行于光栅栅条，然后安装样品，对准成像探测器，拍摄样品在各角度的二维投影像，采集一套样品旋转360度的CT投影数据；

步骤S3，分离出样品旋转180度的吸收投影数据和折射投影数据：

从360度CT投影数据中，选取方向相反的投影光线数据，先扣除无样品的背景光分布，再进行配对组合，分离出样品旋转180度的吸收投影数据和折射投影数据；

若是平行束照明样品，上述配对组合可以直接获得180度的平行束吸收投影数据和折射投影数据；若是扇形束照明样品，还须对扇形束的吸收投影数据和折射投影数据进行重排，经过重排才能获得180度的平行束吸收投影数据和折射投影数据，重排方法和传统CT中扇形束重排方法相同；

步骤S4，根据旋转不变性对180度的平行束折射投影数据进行加权，形成多套满足旋转不变性要求的折射投影数据，然后利用雷当逆变换，重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

所述的光栅剪切相位衬度CT成像数据采集和重建方法，其所述步骤S1中，相干的X射线光束，为平行束或者扇形束；光强随分析光栅相对相位光栅剪切位移而起伏的位移曲线。

所述的光栅剪切相位衬度CT成像数据采集和重建方法，其所述步骤S3中，从360度CT投影数据中，选取方向相反的投影光线数据，先扣除无样品的背景光分布，再进行配对组合，组合方法有两种：①相加，将方向相反的两条光线的投影数据进行相加，分离出样品旋转180度的吸收投影数据；②相减除以相加，分别将方向相反的两条光线的投影数据进行相减和相加，然后将相减组合除以相加组合，分离出样品旋转180度的折射投影数据。

本发明方法根据光栅剪切成像方法和衍射增强成像方法的相似性，将在衍射增强成像方法研究中取得的创新成果推广到光栅剪切成像方法之中，提出适用于人体成像的X射线相位衬度CT数据采集和重建方法。

本发明方法区别于其它X射线相位衬度CT数据采集方法和重建方法的两个显著优势在于，首先它象传统CT数据采集过程一样简单，在采集一套360度投影数据的过程中，除了样品和光源、探测器之间的相对转动之外，无须其它操作；其次它不但可以重建样品的吸收系数和折射率，而且可以重建出衬度更高的折射率一阶导数和折射率二阶导数。因此，本发明提出的相位衬度CT采集方法和重建方法不但简便、易于推广，而且可以重建出衬度更高的折射率导数。

本发明方法不但有别于传统X射线CT数据采集方法和重建方法，而且与目前国际上已经提出的光栅剪切相位衬度CT成像数据采集方法和重建方法不同。与之相比，在投影数据采集过程中，分析光栅只须固定在位移曲线的中值点位置(腰位)，而无须移动，只需要采集一套样品旋转360度的投影数据，投影数据采集可以象传统X射线CT一样迅速；除此之外，本发明提出的重建方法不仅可以重建出样品的吸收系数和折射率，而且可以重建衬度更高的折射率一阶导数和折射率二阶导数。因此，本发明具有将相位衬度三维成像技术推广到医学临床的前景。

附图说明

图1光栅剪切相位衬度CT成像实验装置示意图；

图2投影光线坐标系(r，s，z)和样品坐标系(x，y，z)关系示意图；

图3光栅剪切位移曲线示意图；

图4扇形束坐标系(t，β，z)、投影光线坐标系(r，s，z)和样品坐标系(x，y，z)关系示意图；

图5本发明提出的光栅剪切相位衬度CT成像实验方法流程图。

具体实施方式

本发明基于光栅剪切成像技术，利用分析光栅相对于相位光栅的剪切位移，测量光强随分析光栅剪切位移而变化的位移曲线(形状类似于正弦曲线)，然后分别将分析光栅调整在位移曲线中值点位置(腰位)，并拍摄无样品的背景光分布，将样品转轴与光栅栅条调整平行，采集一套样品旋转360度的CT投影数据。为了分离出样品的吸收投影数据和折射投影数据，从360度的CT投影数据中，选取方向相反的投影光线数据，先扣除无样品的背景光分布，再进行配对组合，相加组合可以获得180度的吸收投影数据，相减组合除以相加组合可以获得180度的折射投影数据。若为平行束照明，组合后的投影数据即为180度的平行束吸收投影数据和折射投影数据。若为扇形束照明，可以通过重排变换，把扇形束的吸收投影数据和折射投影数据变换成180度的平行束吸收投影数据和折射投影数据，重排方法与传统CT中扇形束投影数据的重排方法相同。然后根据旋转不变性要求对180度的平行束折射投影数据其进行加权，形成多套满足旋转不变性要求的180度平行束折射投影数据，然后利用雷当(Radon)逆变换重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

本发明所述的光栅剪切相位衬度CT成像数据采集方法，为了简单有效地获取样品的相位信息，必须将分析光栅调整在位移曲线中值点位置(腰位)，并将样品转轴调整在与光栅栅条平行的位置，采集一套样品旋转360度的CT投影数据。

本发明所述的光栅剪切相位衬度CT成像数据重建方法，必须从360度的CT投影数据中，选取方向相反的投影光线数据进行配对组合和重排变换，形成样品在平行束照射条件下旋转180度的投影数据，再根据旋转不变性对其进行加权，然后利用雷当(Radon)逆变换重建出样品中的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

光栅剪切成像的基本光路示意图如图1所示，成像装置主要由样品转台、相位光栅、分析光栅和成像探测器组成。X射线从样品出射后，在投影光线坐标系(r，s，z)中，其折射角与样品中折射率梯度分量之间的关系为

(1)式和(2)式中的θ_r和θ_s为在(r，s)平面内，分别垂直和平行于光线传播方向的折射角分量；(3)式中的θ_z为(s，z)平面内的折射角分量。(2)式的物理意义为，平行于光线传播方向的折射角分量为零；(2)式在数学上的解释为：因为折射率n在样品外部为1，穿过样品的路径积分

的上限和下限相等，所以积分

必然为零。根据光栅衍射性质，光栅剪切成像只能记录到垂直于栅条方向上的折射角分量，而记录不到平行于栅条方向上的折射角分量，因此本发明只涉及(r，s)平面内的折射角分量θ_r，而不涉及(s，z)平面内的折射角分量θ_z。CT数据采集和重建需要两套坐标系，一套是根据投影光线构造的坐标系(r，s，z)，坐标轴r和坐标轴s分别垂直和平行于投影光线，另一套为固定在样品上的坐标系(x，y，z)，用于描述样品中各种物理量的三维分布；两套坐标系的原点和Z轴重合，它们之间的相对转角为Θ，如图2所示。

本发明基于光栅剪切成像方法和衍射增强成像方法的相似性，根据本专利申请人在衍射增强成像方法研究上取得的创新成果，将衍射增强成像理论推广到光栅剪切成像方法之中，提出光栅剪切成像方程

在(4)式中，I₀为照射在样品上的光强，I为从分析光栅出射的光强，μ为线性吸收系数，θ_r为样品出射光线在(r，s)平面内的折射角分量，D为相位光栅的自成像距离，也是分析光栅和相位光栅之间的距离，x_g为分析光栅和相位光栅之间的剪切位移，R(x_g)为没有样品时，透射光强随分析光栅剪切位移x_g变化的周期性曲线，简称为位移曲线，其中将某一个透射光强最大值的位置选为位移x_g＝0的点，如图3所示。光栅剪切成像方法要求分析光栅的空间周期与相位光栅自成像的空间周期相同，设为p，可得位移曲线的空间周期也为p，相邻两透射光强中值点(两腰位)之间的距离为

也称为位移曲线的半高宽。根据(4)式，当分析光栅位于位移曲线中值点位置(腰位)时，即 $x_g=\frac{p}{4},$ 在样品的折射角分量 ${\mathrm{\&theta;}}_r<\frac{p}{4D}$ 的条件下，光强I与折射角分量θ_r成正比。将(4)式中的 $R(\frac{p}{4}+{\mathrm{\&theta;}}_{r}D)$ 在坐标 $x_g=\frac{p}{4}$ 附近的线性区域展开，得

$R(\frac{p}{4}+{\mathrm{\&theta;}}_{r}D)=R\left(\frac{p}{4}\right)+\frac{\mathrm{dR}\left(\frac{p}{4}\right)}{{\mathrm{dx}}_g}{\mathrm{\&theta;}}_{r}D=R\left(\frac{p}{4}\right)(1+C{\mathrm{\&theta;}}_r),---\left(5\right)$

其中 $C=\frac{D}{R\left(\frac{p}{4}\right)}\frac{\mathrm{dR}\left(\frac{p}{4}\right)}{{\mathrm{dx}}_g}$ 为一个常数。把(1)式和(5)式代入(4)式，得

(6)式和(7)式分别代表在样品旋转时中，方向相反的两根光线穿过样品到达探测器的光强投影数据。从探测器记录的样品旋转360度的CT投影数据中，选取任意两根方向相反的投影光线数据，若进行相加组合，则可以获得180度的吸收投影数据

若将相减组合除以相加组合，则可以获得180度的折射投影数据

若照明样品的光束为平行束，则根据(8)式和(9)式可以分别获得样品旋转180度的平行束吸收投影数据和折射投影数据。若照明样品的光束为扇形束，则需要根据图4中的几何关系，将扇形束重排为平行束。扇束投影可分为等角射线型和等距射线型。等角射线型是指探测器像素等角分布在以光源为中心的圆弧上，等距射线型是指探测器像素等距分布在与中心投影光线垂直的直线上。为简明起见，本专利仅讨论等距射线型情况。在扇形光束照明的情况下，需要引入描述扇形射线的坐标系(t，β，z)，如图4所示。在图4中(a)中，S为光源，绕坐标原点o作半径为l的圆周运动。探测器所在线

也绕o作圆周运动；从射线源发出并通过坐标原点o的中心射线与y轴的夹角为β，探测器

与中心射线垂直，以t表示探测器像素到垂足的距离，角度为β的投影(光强)值是距离t的函数，表示为

(t，β，z)。因为扇形束投影数据和平行束投影数据是等价的完备集，所以经过重排可以将扇形束投影数据变换为平行束投影数据。因此，可以通过扇形射线坐标系(t，β，z)和雷当(Radon)坐标系(r，Θ，z)之间的坐标变换关系，将扇形束投影数据变换为平行束投影数据。根据图4(a)可知，

$\mathrm{\&gamma;}=\arctan \frac{t}{m},---\left(10\right)$

Θ＝β+γ，                                            (11)

$r=\frac{\mathrm{lt}}{\sqrt{m^2+t^2}}.---\left(12\right)$

当光源绕坐标原点o逆时针旋转至图4中(b)位置时，此时光线

正好沿光线

的反方向穿过样品，根据图4(a)中的几何关系可知，

t′＝-t，γ′＝-γ，β′＝β+2γ+π，                  (13)

由(11)、(13)式可得，

Θ′＝β′+γ′＝β+γ+π＝Θ+π。                     (14)

根据(11)式、(12)式和(14)式，可以从360度扇形束投影数据

(t，β，z)获得(8)式和(9)式所表达的180度的平行束吸收投影数据和折射投影数据，有

$\ln \frac{2R\left(\frac{p}{4}\right)I_0}{I(r,\mathrm{\&Theta;},z)+I(r,\mathrm{\&Theta;}+\mathrm{\&pi;},z)}$

$=\ln \frac{2R\left(\frac{p}{4}\right)I_0}{I(\frac{\mathrm{lt}}{\sqrt{m^2+t^2}},\mathrm{\&beta;}+\mathrm{\&gamma;},z)+I(\frac{-\mathrm{lt}}{\sqrt{m^2+t^2}},\mathrm{\&beta;}+\mathrm{\&gamma;}+\mathrm{\&pi;},z)},---\left(15\right)$

$=\ln \frac{2R\left(\frac{p}{4}\right)I_0}{\hat{I}(t,\mathrm{\&beta;},z)+\hat{I}(-t,\mathrm{\&beta;}+\mathrm{\&pi;},z)}$

${\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)=\frac{1}{C}\frac{I(r,\mathrm{\&Theta;},z)-I(r,\mathrm{\&Theta;}+\mathrm{\&pi;},z)}{I(r,\mathrm{\&Theta;},z)+I(r,\mathrm{\&Theta;}+\mathrm{\&pi;},z)}$

$=\frac{1}{C}\frac{I(\frac{\mathrm{lt}}{\sqrt{m^2+t^2}},\mathrm{\&beta;}+\mathrm{\&gamma;},z)-I(\frac{-\mathrm{lt}}{\sqrt{m^2+t^2}},\mathrm{\&beta;}+\mathrm{\&gamma;}+\mathrm{\&pi;},z)}{I(\frac{\mathrm{lt}}{\sqrt{m^2+t^2}},\mathrm{\&beta;}+\mathrm{\&gamma;},z)+I(\frac{-\mathrm{lt}}{\sqrt{m^2+t^2}},\mathrm{\&beta;}+\mathrm{\&gamma;}+\mathrm{\&pi;},z)}.---\left(16\right)$

$=\frac{1}{C}\frac{\hat{I}(t,\mathrm{\&beta;},z)-\hat{I}(-t,\mathrm{\&beta;}+\mathrm{\&pi;},z)}{\hat{I}(t,\mathrm{\&beta;},z)+\hat{I}(-t,\mathrm{\&beta;}+\mathrm{\&pi;},z)}$

$={\widehat{\mathrm{\&theta;}}}_r(t,\mathrm{\&beta;},z)$

因此，在利用扇形束投影数据进行断层重建之前，只须将其重排为平行束投影数据，其余步骤和平行束投影数据重建过程完全相同。

在传统CT中，既可以根据一个物理量的投影重建该物理量，也可以利用根据一个物理量导数(沿垂直于投影光线方向求导)的投影重建该物理量。依据同样的原理，可以利用(8)式中吸收系数的投影重建吸收系数，也可以利用(9)式中折射率导数投影重建折射率。运用傅立叶中心切片定理，可推导出样品吸收系数的重建公式为

$\mathrm{\&mu;}(x,y,z)$

$=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[\ln \frac{2R\left(\frac{p}{4}\right)I_0}{I(r,\mathrm{\&Theta;},z)+I(r,\mathrm{\&Theta;}+\mathrm{\&pi;},z)}*F^{-1}\left(\right|\mathrm{\&rho;}\left|\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr},---\left(17\right)$

参考(1)式，样品折射率的重建公式为

$=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\frac{\left|\mathrm{\&rho;}\right|}{2\mathrm{\&pi;j\&rho;}}\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr}$

下面讨论如何重建折射率导数。CT要求被重建的物理量具有旋转不变性，为此考察(9)式右边积分号内的折射率导数。因为(9)式右边积分号下的折射率导数是在投影光线坐标系中求导的，所以当样品旋转时，折射率导数就会改变数值，即不满足旋转不变性。然而在样品坐标系中求导获得的折射率导数与旋转无关，满足旋转不变性。因此，如果能根据样品坐标系中折射率导数和投影光线坐标系中折射率导数之间的关系，求出两个坐标系中折射率导数投影数据之间的关系，就有可能重建样品坐标系中的折射率导数。根据投影光线坐标系和样品坐标系之间的变换关系

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}\cos \mathrm{\&Theta;}& -\sin \mathrm{\&Theta;}\\ \sin \mathrm{\&Theta;}& \cos \mathrm{\&Theta;}\end{array}\right)\left(\begin{array}{c}r\\ s\end{array}\right),---\left(19\right)$

可得

将(2)式代入(20)式，可得

在(21)式中，θ_x和θ_y分别为在样品坐标系(x，y，z)中获得的折射角分量。

同理有

根据(8)式、(9)式、(21)式和(22)式，可以获得多套满足旋转不变性折射率导数的投影数据，有

运用傅立叶中心切片定理，可以获得样品折射率一阶导数的重建公式为

$\frac{\&PartialD;n(x,y,z)}{\&PartialD;x}=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\cos \mathrm{\&Theta;d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr},---\left(28\right)$

$\frac{\&PartialD;n(x,y,z)}{\&PartialD;y}=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\sin \mathrm{\&Theta;d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr},---\left(29\right)$

样品折射率二维梯度的重建公式为

${\&dtri;}_{\mathrm{xy}}n(x,y,z)=\frac{\&PartialD;n(x,y,z)}{\&PartialD;x}{\stackrel{\&RightArrow;}{e}}_x+\frac{\&PartialD;n(x,y,z)}{\&PartialD;y}{\stackrel{\&RightArrow;}{e}}_y,$

$=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}(\cos \mathrm{\&Theta;}{\stackrel{\&RightArrow;}{e}}_x+\sin \mathrm{\&Theta;}{\stackrel{\&RightArrow;}{e}}_y)\mathrm{d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr}---\left(30\right)$

样品折射率二阶导数的重建公式为

$\frac{{\&PartialD;}^{2}n(x,y,z)}{{\&PartialD;x}^2}=2\mathrm{\&pi;j}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{co}{s}^2\mathrm{\&Theta;d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\mathrm{\&rho;}\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr},---\left(31\right)$

$\frac{{\&PartialD;}^{2}n(x,y,z)}{{\&PartialD;y}^2}=2\mathrm{\&pi;j}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{si}{n}^2\mathrm{\&Theta;d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\mathrm{\&rho;}\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr},---\left(32\right)$

$\frac{{\&PartialD;}^{2}n(x,y,z)}{\&PartialD;x\&PartialD;y}=-2\mathrm{\&pi;j}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\cos \mathrm{\&Theta;}\sin \mathrm{\&Theta;d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\mathrm{\&rho;}\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr};---\left(33\right)$

样品折射率的二维拉普拉斯的重建公式为

${\&dtri;}_{\mathrm{xy}}^{2}n(x,y,z)=\frac{{\&PartialD;}^{2}n(x,y,z)}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^{2}n(x,y,z)}{{\&PartialD;y}^2}$

$=2\mathrm{\&pi;j}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{d\&Theta;}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[{\mathrm{\&theta;}}_r(r,\mathrm{\&Theta;},z)*F^{-1}\left(\right|\mathrm{\&rho;}\left|\mathrm{\&rho;}\right)]\mathrm{\&delta;}(x\cos \mathrm{\&Theta;}+y\sin \mathrm{\&Theta;}-r)\mathrm{dr}.---\left(34\right)$

根据图1所示的光栅剪切成像基本光路，在光路上依次排列样品转台、相位光栅、分析光栅和成像探测器；沿着投影光束方向将分析光栅调整在相位光栅自成像的距离，再沿着垂直于投影光束的方向移动分析光栅，测量光强随分析光栅剪切移动而变化的位移曲线，然后将分析光栅固定在位移曲线中值点的位置(腰位)上，并拍摄无样品的背景光分布，调整样品转轴方向，使其与光栅栅条平行，安装样品，对准成像探测器，对样品进行360°旋转扫描，拍摄样品在各角度的二维投影像，采集一套样品旋转360°的投影数据。

从360度的CT投影数据中，选取方向相反的投影光线数据，先扣除无样品的背景光分布，再进行配对组合，将吸收投影信号和折射投影信号分离，经过重排变换，形成样品旋转180度的平行束吸收投影数据和折射投影数据，再根据旋转不变性对其进行加权，形成多套满足旋转不变性要求的平行束180度投影数据，然后根据(17)式、(18)式、(28)式、(29)式、(31)式、(32)式和(33)式，重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

图5是光栅剪切相位衬度CT成像实验方法流程图。其具体步骤如下：

步骤S1，准备光栅剪切成像条件：相干的X射线光束(平行束或者扇形束)照射相位光栅，根据光学中的泰保(Talbot)效应，在相位光栅后面会产生自成像，自成像的空间频率是相位光栅的空间频率的一倍或者二倍，在自成像位置，放置空间频率和自成像空间频率相同的分析光栅(又称为吸收光栅)，沿着垂直于光束传播的方向，移动分析光栅，同时利用光强探测器测得光强随着分析光栅位置移动而起伏的位移曲线(位移曲线形状类似于正弦曲线)，并拍摄无样品的背景光分布。

步骤S2，采集样品旋转360度的折射衬度CT投影数据：将分析光栅位置调整在位移曲线的中值点，在相位光栅前安装样品转台，调整样品转台旋转轴方向，使其平行于光栅栅条，然后安装样品，对准成像探测器，拍摄样品在各角度的二维投影像，采集一套样品旋转360度的CT投影数据。

步骤S3，分离出吸收投影数据和折射投影数据：从360度CT投影数据中，选取方向相反的投影光线数据，先扣除无样品的背景光分布，再进行配对组合，分离出180度的吸收投影数据和折射投影数据。组合方法有两种：①相加，②相减，③相减除以相加。将方向相反的两条光线的投影数据进行相加，可以分离出样品吸收投影数据。分别将方向相反的两条光线的投影数据进行相减和相加，然后将相减组合除以相加组合，可以分离出样品折射投影数据。若是平行束照明样品，上述配对组合可以直接获得180度的平行束吸收投影数据和折射投影数据。若是扇形束照明样品，由于扇形束投影数据和平行束投影数据对于断层重建而言是等价的完备集，所以通过重排变换，可以把扇形束投影数据变换成平行束投影数据。为此需要对扇形束的吸收投影数据和折射投影数据进行重排，重排出180度的平行束吸收投影数据和折射投影数据，重排方法和传统CT中扇形束重排方法相同。

步骤S4，根据旋转不变性对折射投影数据进行加权，形成多套满足旋转不变性要求的折射投影数据，然后利用雷当(Radon)逆变换，重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。因为吸收系数和折射率都是标量，满足旋转不变性条件，所以不需要加权。折射率一阶导数不是标量，需要分别为

和

的投影数据乘以权重因子cosΘ和sinΘ，其中(x，y，z)固定在样品上的坐标系，样品围绕z轴旋转，Θ为样品转角。同理折射率二阶导数也不是标量，需要分别为

和

的投影数据乘以权重因子cos²Θ、sin²Θ和cosΘsinΘ。根据(17)式、(18)式、(28)式、(29)式、(31)式、(32)式和(33)式，重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

Claims (3)

1.一种光栅剪切相位衬度CT成像数据采集和重建方法，其特征在于：在分析光栅处于位移曲线中值点位置，样品转轴与光栅栅条平行的条件下，采集一套样品旋转360度的CT投影数据；从360度的CT投影数据中选取方向相反的投影光线数据进行配对组合，分离出样品旋转180度的吸收投影数据和折射投影数据，再经重排变换，形成样品旋转180度的平行束吸收投影数据和折射投影数据；根据旋转不变性要求对折射投影数据进行加权，获得多套满足旋转不变性要求的180度平行束投影数据，再根据雷当逆变换重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布；

具体步骤如下：

步骤S1，准备光栅剪切成像条件：

用相干的X射线光束照射相位光栅，根据光学中的泰保效应，在相位光栅后面会产生自成像，自成像的空间频率是相位光栅的空间频率的一倍或者二倍，在自成像位置，放置空间频率和自成像空间频率相同的分析光栅，沿着垂直于光束和光栅栅条的方向，移动分析光栅，同时利用光强探测器测得光强随着分析光栅位置移动而起伏的位移曲线，并拍摄无样品的背景光分布；

步骤S2，采集样品旋转360度的相位衬度CT投影数据：

将分析光栅位置调整在位移曲线的中值点，在相位光栅前安装样品转台，调整样品转台旋转轴方向，使其平行于光栅栅条，然后安装样品，对准成像探测器，拍摄样品在各角度的二维投影像，采集一套样品旋转360度的CT投影数据；

步骤S3，分离出样品旋转180度的吸收投影数据和折射投影数据：

从360度CT投影数据中，选取方向相反的投影光线数据，先扣除无样品的背景光分布，再进行配对组合，分离出样品旋转180度的吸收投影数据和折射投影数据；

若是平行束照明样品，上述配对组合直接获得180度的平行束吸收投影数据和折射投影数据；若是扇形束照明样品，还须对扇形束的吸收投影数据和折射投影数据进行重排，经过重排才能获得180度的平行束吸收投影数据和折射投影数据，重排方法和传统CT中扇形束重排方法相同；

步骤S4，根据旋转不变性对180度的平行束折射投影数据进行加权，形成多套满足旋转不变性要求的折射投影数据，然后利用雷当逆变换，重建出样品的吸收系数、折射率、折射率一阶导数和折射率二阶导数的三维分布。

2.如权利要求1所述的光栅剪切相位衬度CT成像数据采集和重建方法，其特征在于：所述步骤S1中，相干的X射线光束，为平行束或者扇形束；光强随分析光栅相对相位光栅剪切位移而起伏的位移曲线。

3.如权利要求1所述的光栅剪切相位衬度CT成像数据采集和重建方法，其特征在于：所述步骤S3中，从360度CT投影数据中，选取方向相反的投影光线数据进行配对组合，组合方法有两种：①相加，将方向相反的两条光线的投影数据进行相加，分离出样品旋转180度的吸收投影数据；②相减除以相加，分别将方向相反的两条光线的投影数据进行相减和相加，然后将相减组合除以相加组合，分离出样品旋转180度的折射投影数据。

Images