CN103098095B - 微分相位对比成像中的规则化的相位恢复 - Google Patents

Abstract

本发明涉及对象（108）的微分相位对比成像。当从微分相位对比图像数据重建图像信息时，可能产生条纹状伪影（502）。所述伪影（502）可能实质地降低重建图像数据的易读性。因此，去除或者至少抑制所述伪影（502）可能是有益的。因此，提供了用于相位对比成像中的规则化的相位恢复的方法（400），包括：接收（402）对象（108）的微分相位对比图像数据；生成（404）对象（108）的重建图像数据；并呈现（406）对象（108）的重建图像数据。微分相位对比图像和重建图像数据包括具有第一维度和第二维度的二维图像数据结构。生成重建图像数据包括在数据结构的第一维度和第二维度之一中对图像数据的积分，在所述数据结构的第一维度和第二维度的另一个中确定梯度算子，并且采用数据结构以重建图像数据实现重建图像数据内的伪影（500）的降低。

Description

微分相位对比成像中的规则化的相位恢复

技术领域

本发明总体涉及X射线成像技术。

更具体而言，本发明涉及微分相位对比成像。

具体而言，本发明涉及一种用于相位对比成像中的规则化（regularized）的相位恢复（retrieval）的方法、一种用于相位对比成像中的规则化的相位恢复的设备、一种包括根据本发明的设备的X射线系统、一种根据本发明的设备在X射线系统和CT系统之一中的用途、一种计算机可读介质以及一种程序单元。

背景技术

在透射X射线图像采集中，待检查的对象，例如患者，被布置在例如X射线管的X射线生成装置与X射线探测器之间。

源自X射线生成装置的X射线辐射贯穿待检查的对象并且随后到达用于采集图像信息的X射线探测器，所采集的图像信息之后被重建为X射线图像以供呈现。对象的内部结构，例如，组织结构，提供对贯穿所述对象的X辐射的空间衰减。相应地，X射线探测器配准空间衰减的X射线辐射。

某些对象可能仅对X射线辐射衰减较小的程度，或者可能相当均匀地衰减X射线束，造成所采集的X射线图像具有低对比度。

然而，即使对象仅对贯穿对象的X射线束有少量的衰减，同一对象对X射线辐射的波前的相位的影响的程度也可能相当大。

相应地，可以采用相位对比成像以获得对象的具有增强的对比度的图像信息。

在相位对比成像中，X射线源连同邻近所述X射线源布置的所谓的源格栅元件生成至少部分空间相干的X射线辐射。贯穿对象的相干的X射线允许相位信息的之后的恢复。

由于波的相位可能不是直接测量的，采用另一格栅元件，即所谓的相位格栅，其被布置在待检查的对象与X射线探测器之间。相位格栅允许通过多个波的干涉将相移转换为强度调制，所述强度调制然后能够由X射线探测器探测。

然而，通过仅采用相位格栅生成的干涉模式（pattern）对于当前的X射线探测器而言可能过小而不能够被精确地探测，因为X射线探测器的空间分辨率不足。在本文中，可以采用额外的格栅元件，即所谓的分析器格栅，其被邻近X射线探测器布置在相位格栅元件与X射线探测器之间。分析器格栅提供干涉模式，其足够大从而能够被当前X射线探测器探测。

为了获得适当的相位对比图像信息，执行相位步进以获得多个相位对比投影。在相位步进中，源格栅元件、相位格栅元件和分析器格栅元件之一被相对于其他格栅和X射线探测器元件横向位移相应格栅跨度（pitch）的分数，例如，例如相位格栅的格栅跨度的四分之一、六分之一、八分之一。将图像采集和横向位移重复例如四次、六次或八次，以采集多个相位对比投影，共同构成相位步进间隔。

在微分相位对比成像中，X射线探测器探测并由此测量垂直于格栅元件的格栅方向的，即，垂直于屏蔽区域和格栅结构的沟槽区域的延长线的一阶导数。由于格栅结构，所采集的图像信息可以被认为是高度对称的，其中边缘在一方向上，即，在垂直于格栅方向的方向上，尤其被增强。平行于所述格栅方向的方向可能未被增强。

为了重建图像信息，沿微分的线，即，垂直于格栅结构的积分流程，可以得到具有条纹状伪影或者由于噪声或其他误差的条纹的图像信息，其局部地布置在微分图像信息中并且其可以沿图像数据的记忆的线传播。

由此，可能有益的是提供用于降低甚至移除所述条纹状伪影的装置。

在Weitkamp T.、Diaz A.、David C.等人的“X-ray phase imaging with agrating interferometer”，Optics Express6296,8.August2005/vol.13,no.16，以及Bartl P.、Durst J.、Haas W.等人的“Simulation of X-ray phase-contrast computedtomography of a medical phantom comprising particle and wave contributions”Proc of SPIE vol.762276220Q-1这两篇文献中都描述了相位对比成像。

在Numer.Math.69:25-31（1994）上的Engl H.W.和Grever W.的“Using the L-curve for determining optimal regularization parameters”一文中描述了确定规则化参数。

ZHIHUA QI ET AL在PROCEEDINGS OF SPIE，第7258卷，2009年1月1日（2009-01-01），第72584A-1-72584A-8页，XP55013373，ISSN:0277-786X，DOI:10.1117/12.813861上的“A novel method to reduce data acquisition time in differential phasecontrast:computed tomography using compressed sensing”一文描述了基于采用在两个方向的稀疏变换的压缩感测架构的重建方法。

发明内容

本发明的一方面可以被视为采用已知的规则化方法来解决不良姿态（ill-posed）问题以降低图像伪影。

相应地，提供了根据独立权利要求的一种用于相位对比成像中的规则化的相位恢复的方法、一种用于相位对比成像中的规则化的相位恢复的设备，一种包括根据本发明的设备的X射线系统、根据本发明的设备在X射线系统和CT系统之一中的使用、一种计算机可读介质以及一种程序单元。

从属权利要求能够得到优选实施例。

本发明的这些和其他方面能够参考上下文中描述的实施例变得显而易见并得以阐述。

将参考以下附图描述本发明的示范性实施例。

附图中的图示说明是示意性的。在不同的附图中，以相似或相同的参考标记提供相似或相同的元件。

附图未按照比例绘制，然而可以描绘定性性质。

附图说明

图1示出了根据本发明的X射线系统的示范性实施例；

图2示出了根据本发明的微分相位对比成像系统的示范性实施例；

图3示出了根据本发明的用于相位对比成像中的规则化的相位恢复的设备的示范性实施例；

图4示出了根据本发明的用于相位对比成像中的规则化的相位恢复的方法的示范性实施例；

图5a-c示出了根据本发明的采用微分相位对比图像数据的图像信息的示范性重建。

参考标记列表

100 X射线系统

102 X射线生成装置/X射线源

104 X射线探测器

106 机架

108 对象

110 支撑物

112 处理系统

114 X射线辐射

200 微分相位对比成像系统

202 源格栅

204 相位格栅G1

206 分析器格栅G2

208 致动器元件

210a、b 波前

300 微处理器

302 存储元件

304 显示元件

400 用于规则化的相位恢复和相位对比成像的方法

402 接收微分相位对比图像数据

404 生成重建的图像数据

406 呈现重建的图像数据

500 条纹状伪影

具体实施方式

本发明应用用于解决不良姿态问题的规则化的方法，以解决当从采用微分相位对比成像系统采集的微分相位对比图像数据来重建X射线图像时的伪影抑制问题。

对于以下的解释说明，不失一般性地，假设所采集的二维图像信息是包括NxN矩阵图像结构的图像信息。换言之，二维图像包括图像行与图像列同样多的像素，由此得到二次方图像。因此，用于采集所述图像数据而采用的X射线探测器元件可以例如包括512x512或者1024x1024的单个个体像素元件。

此外，假设待重建的图像信息包括Nx(N+1)的矩阵形式。

然而，应当理解，以下解释还可以应用于具有十分不同的像素结构的图像数据，并且本领域技术人员可以容易地调整以下教导以准确地实现图像数据大小或形状。

对于以下解释说明，我们假设条纹状伪影作为源自相对于图像信息竖直对准的格栅结构的水平条纹。此外，不失一般性的，假设微分的方向是沿着图像的行。图像数据的行由i表示，而列由j表示。

测量到的微分图像数据被称为二维图像数据d_i,j，而期望的二维图像信息，即重建的二维图像信息，被表示为g_i,j。

对于规则化方法的应用而言，逐行积分问题被再次变换构成线性重建问题。微分操作可以被视为前向（forward）问题，其可以根据方程1被表示为g与d之间的关系。

d_i,j=g_i,j+1-g_i,j 方程1

根据方程1，给出了g与d之间的线性关系。因此，二维图像数据g和d可以被表示为一维矢量，而前向操作可以根据方程2被用方程表示为矩阵A。

d=A·g 方程2

由于图像信息g是要从测量的图像数据d或测量的矢量d重建的图像信息，重建操作可以被视为从所测量的矢量d来估计g。

相对于矢量记法，针对d和g的单个脚本指示，例如d_k，可以被用于识别相应矢量d和g的单个元件。

所述指示可以根据关系k(i,j)=Ni+j或者反之j(k)=k mod N，i(k)=(k/N)的整数部分来映射。那么矩阵A是如下稀疏矩阵：

为了降低条纹状伪影，例如，水平条纹伪影，可以采用梯度算子，其处罚，即，降低与二维图像数据的积分的维度不同，即，与执行积分的维度不同的第二维度中的梯度。

相应地，可以根据如下方程定义梯度算子：

${\left({\▿}_{z}g\right)}_{i,j}=g_{i+1,j}-g_{i,j}$ 方程3

因此，在二维图像数据的探测器行i之内执行积分操作的情况下，梯度算子对二维图像数据的每列j进行操作。

如前向算子矩阵A，梯度算子可以被理解为是对一维矢量g进行操作的线性算子。期望的图像信息g由此可以通过根据方程4对最小化问题进行求解来进行重建。

${\left|\right|A\·g-d\left|\right|}^2+\mathrm{\λ}\·{\left|\right|{\▿}_z\·g\left|\right|}^2=\min$ 方程4

规则化参数λ可以被用于控制应当获得的平滑度，即，条纹伪影降低的量。

可能需要通过正在视觉观察根据规则化参数λ的重建图像数据的用户对规则化参数λ的内部调整来经验性地找到λ。所述算子因此可以在内部控制规则化参数直到获得了优选的重建图像。规则化参数λ的可能的值是实正值。λ还可以通过采用L曲线或者通过在EnglH.W描述的Morozov的差别（discrepancy）原理来确定。

方程4可以被视为惩罚最大似然重建算法。为了获得更为一般化的形式，可以采用两种一般化。

首先，所测量的数据d_k可以被视为高斯分布随机变量，其中，方差为σ_k ²。甚至更为一般化的，噪声协方差矩阵C_kl可以与对应于σ_k ²的对角要素C_kk一起采用。

第二，z方向的梯度可以由势函数来不同地进行惩罚。

为了获得最小化问题的唯一解，所采用的势函数可能需要是凸状的。

对应的势函数可以是二次惩罚、根据方程5的Huber惩罚或Huber损失函数或者根据方程6的一般化的高斯Markov随机场之一：

方程5

${\mathrm{\ψ}}_{p,q,c}\left(x\right)=\frac{{\left|x\right|}^p}{1+{\left|\frac{x}{c}\right|}^{p-q}}$ 方程6

其中，Huber惩罚中的σ>0以及一般化高斯Markov随机场中的p、q和c（其中，p≥q≥1和c>0）是经验参数，其是在前面选择的。

通过采用对应的一般化，规则化的恢复的相位前可以通过根据方程7对最小化问题求解来获得。

${(A\·g-d)}^T\·C^{-1}\·(A\·g-d)+\mathrm{\λ}\·\underset{k}{\mathrm{\Σ}}{w}_k\mathrm{\ψ}\left({({\▿}_z\·g)}_k\right)=\min \·$ 方程7

在方程7中，额外的权重因子w_k可以局部地强调或去强调由于惩罚在z方向的平滑，并且是实非负值。

方程4以及方程7两者可以被视为试图针对每个探测器行重建积分的额外常数，因为N+1值是从仅N个数据样本针对每行估计的，所述常数可以被考虑为是条纹状伪影的原因。

在待检查的对象被整体匹配到X射线辐射的单个扇形中的情况下，矩阵g_i,j还可以是NxN大小矩阵，其中积分常量可以被固定为0。具体而言，为了获得微分数据d，所测量的数据是指所谓的空白扫描，即，无对象的扫描。通过将具有对象的扫描参考空白扫描，所述处理变得仅对对象影响的波前的改变敏感，而不是对波前自身的相位敏感。如果对象完全匹配到扇形中，每个探测器行中的第一像素已知要受到对象的影响，由此，与空白扫描相比相位的改变为零。

惩罚，即，降低在例如z方向的方向上的梯度可以被考虑为对于抑制特别是横向伪影的条纹状伪影而言非常有效。条纹状伪影可以被考虑为创建在整幅图像上的图像，其中惩罚快速加和到实质的值。相应地，通过仅调节单个估计的测量值可以去除条纹伪影，即，数据惩罚项可能不与惩罚项一样快地改变。这能够被容易地认为这种情况的范例，其中，数据d中的噪声由离群值（outlier）支配，即，由数据d₂₁中的异常大的噪声值支配。通过标准的相位积分，即，通过根据来估计g，通过横贯整个行的和传播单个外轮廓，并且由此这一外轮廓将针对大数量的探测器像素在z中的梯度做出贡献。

应当理解，从反转的场获知的方法的应用，不良姿态问题既不是已知的也不是降低图像信息内的伪影的证据，因为实际重建是算术的，不是不良姿态问题，而是良好姿态的问题。

现在参考图1，描绘了根据本发明的X射线系统的示范性实施例。

X射线系统100，示范性描绘为CT系统100，包括X射线生成装置102，例如X射线管；和X射线探测器104，其包括个体探测器像素元件116。X射线探测器元件104被示范性地描绘为探测器像素元件116的二维阵列。

X射线生成装置102和X射线探测器104两者被彼此相对地布置在机架106上。机架106允许X射线生成装置102和X射线探测器102关于一轴进行旋转。X射线辐射114是由X射线生成装置102生成的并其随后到达X射线探测器104。X射线生成装置102还可以包括源格栅元件G₀202，其在图1中未描绘，用于生成至少部分空间相干的X射线辐射114，以用于微分相位对比成像。像例如相位格栅元件G₁204的其他格栅元件以及分析器格栅元件G₂206在图1中未描绘。

待检查的对象108被布置在支撑物110上，所述支撑物适于允许将对象108移动和放置在X射线束114之内。

处理系统112通信地耦合到X射线系统100，用于控制图像采集流程，用于重建图像数据和/或用于重建的图像信息的随后的呈现，例如显示。

现在参考图2，描绘了根据本发明的微分相位对比成像系统的示范性实施例。

源格栅元件G₀202被示意性地布置在X射线生成装置102的附近，用于生成至少部分空间相干的X射线辐射114。源格栅元件G₀202被布放置成与相位格栅元件G₁204空间分离距离l。

对象108被布置在X射线生成装置102与相位格栅元件G₁204之间。设备200，相位对比成像系统200，还包括与相位格栅元件G₁204空间分离距离D的分析器格栅元件G₂206。分析器格栅元件G₂206被布置在包括个体探测器像素元件116的X射线探测器元件104的附近。

波前210a到达对象108处，具有根据相干X射线辐射的均匀相位关系。在通过对象108之后，均匀相位关系可能已经受到影响，如由波前210b所描绘的。在通过相位格栅元件G₁204和分析器格栅G₂206之后，干涉模式被投影在X射线探测器元件104及其个体探测器像素元件116上。

致动器元件208被示范性地布置在分析器格栅元件G₂206处，用于分析器格栅元件G₂206关于相位格栅元件G₁204的跨度p的分数的位移，用于在相位步进期间采集个体相位对比投影，即，在单个相位步进采集间隔中具有个体相位步进状态。

致动器元件208还可以被布置在相位格栅元件G₁204和源格栅元件G₀202之一处。

现在参考图3，描绘了根据本发明的用于相位对比成像中的规则化的相位恢复的设备的示范性实施例。

处理系统112示范性地包括处理元件300，其被附接到存储元件302和显示元件304。

处理系统112通信地与X射线系统100耦合，用于控制X射线系统，用于重建所采集的图像信息和/或用于呈现所重建的图像信息，例如，在显示器304上显示重建的图像信息。

然而，呈现所重建的图像信息还可以被理解为向其他处理系统或显示系统提供所述重建的图像信息，和/或本地存储重建的图像信息，例如，用于长期存储和存档。

处理元件300适于执行用于相位对比成像中的规则化的相位恢复的方法400，具体而言，用于根据本发明的方法400重建所采集的或者所接收的对象的相位对比图像数据。

现在参考图4，描绘了用于根据本发明的在相位对比成像中的规则化的相位恢复的方法的示范性实施例。

图4示出了用于相位对比成像中的规则化的相位恢复的方法400，包括接收402对象的相位对比图像数据，生成404对象的重建的图像数据，并呈现对象的重建的图像数据。

重建的图像数据的生成可以包括采用根据本发明的规则化方法，用于重新组合图像数据，具体而言，用于去除或抑制像例如重建图像数据中的条纹状伪影的伪影。

呈现对象108的重建图像数据由此可以是向X射线系统100的用户显示图像信息。然而，呈现还可以包括存储重建的图像数据，以供重建图像数据的之后的显示、归档，由此以供图像数据的永久或暂时存储，并且所述呈现还可以包括向用于例如显示、存储或进一步处理的其他处理系统来提供图像数据。

现在参考图5a-c，描绘了根据本发明的采用微分相位对比图像数据的图像信息的示范性重建。

在图5a中，呈现了二维微分相位对比图像数据。图5a的图像数据描绘了如在微分相位对比成像中由X射线探测器104采集的图像的个体像素的微分图像关系。而个体列是指索引j，个体行是指索引i。

在图5a中，个体像素的灰度值充分描绘了图像数据相对于邻近像素的差异。例如，在图5a的环状对象108的背景中，呈现了均匀灰度值，例证了图像信息相对于邻近像素基本没有改变。相应地，在图5b和5c中，对象108的对应背景被描绘为黑色。对象108的左半球处的亮像素值描绘了图像信息的相当陡峭的改变，或者正的或者负的。对象108的右半球处的暗像素值例证了在其他方向的基本相似陡峭的改变，由此或者为负的或者为正的。图5b和5c中的图像信息由此能够通过积分获得，例如，沿图5a的个体行（i=常量）的积分。

图5b描绘了采用通用重建算法，例如，积分运算，根据图5a的微分相位对比成像信息的重建图像。条纹状伪影500，在图5中呈现示例性描绘的横向条纹。图5b的图像信息是通过图5a的图像信息的逐行积分从图5a的图像信息获得的。

现在，当采用根据本发明的方法400时，在重建图像信息时，可以获得根据图5c的X射线图像。对于图5c的图像信息，除了逐行积分之外，采用在z方向，即，相对于列（j=常量）的梯度算子惩罚梯度。通过采用根据本发明的梯度算子，从图5c中呈现的重建图像数据基本去除了条纹状伪影500。因此，通过采用根据本发明的梯度算子，可以被降低或者甚至去除条纹状伪影500。

应当注意到，术语“包括”不排除其他元件或步骤，并且“一”或“一个”不排除多个。同样地，与不同实施例相关联描述的元件可以进行组合。

应当注意到，权利要求中的附图标记不应当被解释为限制权利要求的范围。

Claims (15)

1.一种用于相位对比成像中的规则化的相位恢复的方法(400)，包括：

接收(402)对象(108)的微分相位对比图像数据d_i,j；

生成(404)所述对象(108)的重建图像数据g_i,j；以及

呈现(406)所述对象(108)的重建图像数据g_i,j；

其中，所述微分相位对比图像数据和所述重建图像数据包括具有第一维度i和第二维度j的二维数据结构；

其中，生成重建图像数据包括在所述数据结构的所述第一维度i和所述第二维度j之一中对图像数据的积分；

其中，在所述数据结构的所述第一维度i和所述第二维度j的另一个中确定梯度算子并且

其中，采用所述梯度算子以用于重建图像数据。

2.根据权利要求1所述的方法，

其中，所述重建图像数据是通过采用如下方程生成的：

d_i,j＝g_i,j+1-g_i,j。

3.根据权利要求1或2所述的方法，

其中，所述梯度算子是由如下方程确定的：

${\left({\▿}_{Z}g\right)}_{i,j}=g_{i+1,j}-g_{i,j}.$

4.根据权利要求1或2所述的方法，

其中，所述重建图像数据是通过对如下方程求解来重建的：

$\left|\right|A\·g-d|{|}^2+\mathrm{\λ}\·||{\▿}_Z\·g|{|}^2=min,$

其中，A是前向算子矩阵，并且λ是规则化参数。

5.根据权利要求1或2所述的方法，

其中，所述重建图像数据是通过对如下方程求解来重建的：

${(A\·g-d)}^T\·C^{-1}\·(A\·g-d)+\mathrm{\λ}\·\underset{k}{\Σ}{w}_k\mathrm{\ψ}\left({\left({\▿}_Z\·g\right)}_k\right)=min,$

其中，A是前向算子矩阵，λ是规则化参数，ψ是势函数，w_k是权重因子，C是噪声协方差矩阵，并且k(i，j)＝Ni+j。

6.根据权利要求5所述的方法，

其中，ψ是二次惩罚、Huber惩罚以及一般化高斯Markov随机场的一个势函数。

7.一种用于相位对比成像中的规则化的相位恢复的设备(112)，包括：

存储元件(302)，其用于存储接收的对象的微分相位对比图像数据d_i,j；

处理元件(300)，其用于生成所述对象的重建图像数据g_i,j；

其中，所述设备(112)适于执行根据前述权利要求之一所述的方法(400)。

8.一种包括根据权利要求7所述的设备(112)的X射线系统(100)。

9.一种根据权利要求7所述的设备(112)在X射线系统(100)中的用途。

10.一种用于相位对比成像中的规则化的相位恢复的装置，包括：

用于接收(402)对象(108)的微分相位对比图像数据d_i,j的模块；

用于生成(404)所述对象(108)的重建图像数据g_i,j的模块；以及

用于呈现(406)所述对象(108)的重建图像数据g_i,j的模块；

其中，所述微分相位对比图像数据和所述重建图像数据包括具有第一维度i和第二维度j的二维数据结构；

其中，生成重建图像数据包括在所述数据结构的所述第一维度i和所述第二维度j之一中对图像数据的积分；

其中，在所述数据结构的所述第一维度i和所述第二维度j的另一个中确定梯度算子并且

其中，采用所述梯度算子以用于重建图像数据。

11.根据权利要求10所述的装置，

其中，所述重建图像数据是通过采用如下方程生成的：

d_i,j＝g_i,j+1-g_i,j。

12.根据权利要求10或11所述的装置，

其中，所述梯度算子是由如下方程确定的：

${\left({\▿}_{Z}g\right)}_{i,j}=g_{i+1,j}-g_{i,j}.$

13.根据权利要求10或11所述的装置，

其中，所述重建图像数据是通过对如下方程求解来重建的：

$\left|\right|A\·g-d|{|}^2+\mathrm{\λ}\·||{\▿}_Z\·g|{|}^2=min,$

其中，A是前向算子矩阵，并且λ是规则化参数。

14.根据权利要求10或11所述的装置，

其中，所述重建图像数据是通过对如下方程求解来重建的：

${(A\·g-d)}^T\·C^{-1}\·(A\·g-d)+\mathrm{\λ}\·\underset{k}{\Σ}{w}_k\mathrm{\ψ}\left({\left({\▿}_Z\·g\right)}_k\right)=min,$

其中，A是前向算子矩阵，λ是规则化参数，ψ是势函数，w_k是权重因子，C是噪声协方差矩阵，并且k(i，j)＝Ni+j。

15.根据权利要求14所述的装置，

其中，ψ是二次惩罚、Huber惩罚以及一般化高斯Markov随机场的一个势函数。