论文阅读：A digitally reconstructed radiograph algorithm calculated from first principles

最新推荐文章于 2024-05-26 19:58:52 发布

A2333fun

最新推荐文章于 2024-05-26 19:58:52 发布

阅读量946

点赞数 1

分类专栏：论文阅读文章标签：算法深度学习机器学习

本文链接：https://blog.csdn.net/qq_42042757/article/details/120012693

版权

论文阅读专栏收录该内容

21 篇文章 8 订阅

订阅专栏

A digitally reconstructed radiograph algorithm calculated from first principles

基于第一性原理计算的数字重建射线照相算法

最近在研究DRR的生成方式，有兴趣一起讨论加q734405202

**目的：**开发一种算法，用于计算符合真实锥形束CT (CBCT)投影的真实数字重建x线照片(drr)，无需人工调整。

CBCT就是Cone beam CT的简称，即锥形束CT。顾名思义是锥形束投照计算机重组断层影像设备，其原理是X线发生器以较低的射线量（通常球管电流在10毫安左右）围绕投照体做环形DR（数字式投照）。然后将围绕投照体多次（180次-360次，依产品不同而异）数字投照后“交集”中所获得的数据在计算机中“重组，reconstruction”后进而获得三维图像。CBCT获取数据的投照原理和传统扇形扫描CT是完全不同的，而后期计算机重组的算法原理有类似之处。

**方法：**作者利用不同材料的CBCT投影x线片去测量的衰减数据，得到一个函数，这个函数可以将将CT数转换为线性衰减系数(LAC)。首先从衰减数据中去除散射、光束硬化和遮蔽眩光的影响。利用这个转换函数，作者通过CT沿着连接辐射源和探测器像素点的射线计算LAC的线积分，使用射线追踪算法，生成原始drr。通过后处理，散射、光束硬化和遮蔽眩光的影响被包括在drr中。

结果：作者将实际CBCT投影与经过所有校正(散射、光束硬化和遮蔽眩光)产生的drr以及未校正的drr进行了比较。通过投影与drr的视觉比较、像素强度比较、强度直方图比较以及drr与投影像素强度的相关图来评估算法的准确性。总的来说，完全修正的算法在精度上比未修正的算法有一个很小但不是微不足道的改进。作者还研究了基于测量和基于计算的方法来确定梁硬化校正，并发现基于计算的方法更优越，因为它可以解释不均匀的蝶形滤波器厚度。作者对该算法的速度进行了基准测试，发现在0.5 mm的射线步长下，对于全探测器和CT分辨率，该算法产生drr的时间约为0.35秒。

**总结：**作者演示了一种DRR算法，该算法是根据散射、光束硬化和遮蔽眩光的基本原理计算的，以产生准确的DRR。该算法计算效率高，是要求基于drr和投影匹配的数据保真度项的迭代CT重建技术的良好候选算法。

1、介绍

我们正在开发一种前向迭代投影匹配(FIPM)算法，用于原始锥形束计算机断层成像(CBCT)投影和扇形束计算机断层成像(FBCT)图像的可变形2D/3D配准。【1 - 4】这个过程需要从FBCT计算数字重建x线照片(drr)，以与实际的锥形束x线投影进行比较。解剖特征的配准，以及drr和投影比较的定量方式(即匹配度量，或价值函数)，决定了两个图像匹配的准确性。

例如，放疗应用中的刚性2D/3D配准通常聚焦于高对比度的骨骼特征，而可变形配准算法必须对低对比度的软组织特征敏感。在drr到投影匹配的其他应用中，我们注意到同时代数重建技术和先验图像约束遥感【5,6】。每一个FIPM应用程序对DRR质量都有自己的要求。

【7-9】通过投影匹配到DRRs的刚性2D/3D配准问题有着悠久的研究历史，在此期间，人们研究了许多不同的图像匹配度量。虽然我们的兴趣是可变形2D/3D配准，但从刚性配准中获得的见解可以为我们提供有关此类优化参数的有用指导原则。
【4,10,11】文献中记录的长期经验表明，每个指标都有优缺点。平方强度差之和（SSD）具有对噪声鲁棒性的优点，并且对于梯度驱动的迭代搜索过程容易形成，但是需要DRR和投影中绝对图像强度的密切对应。【10】互信息的许多公式允许在比较图像强度中有更大的可变性，但其有效性的大多数研究都涉及3D图像。由于与3D/3D配准相比，可用数据样本的数量相对较少，已经显示出互信息对于2D图像匹配的健壮性较差。当有噪音时，这一弱点会变得更加严重。它的捕获范围更小，界限更清晰，这使得它在梯度驱动的迭代搜索中不那么健壮。互信息迭代代价函数的梯度的确定也比SSD算法更加困难。

在可变形FIPM的开发中，我们首选梯度搜索，因此选择使用差异的平方和作为我们的图像匹配度量。我们已经通过仿真表明，这个公式对drr和投影图像之间的对比度不匹配很敏感，但对噪声相对不敏感。这就是目前工作的动机。**我们的目标是仅从图像形成的基本物理原理计算drr，以准确复制CBCT投影中的绝对图像强度;即，从第一原则计算drr。**我们避免使用强度归一化或任何其他人为的数值强度调整，这可能会使配准倾向于不准确的变形。我们的目的不是最大化drr中的对比度，而是再现投影中的内容。噪声对FIPM的影响很小，而且相对容易合并，因此不在这里考虑。

许多研究人员已经描述了涉及DRR生成过程各个方面的算法。特别是x射线体绘制组件，在几何真实感、空间分辨率和计算效率方面受到了很大的关注。【12】西登介绍了一种广泛使用的算法，可以有效地计算出通过计算机断层扫描(CT)体积的精确放射路径。【13】Russakoffet开发了“衰减场”的概念，用于快速DRR渲染。【14】Moore使用铅笔光束而不是射线作为投影仪，以增加DRR真实感和避免混叠效果。【15】Birkfellneret使用了一种“摆动飞溅”技术来快速生成DRR。【16-18】其他技术包括剪切变形体绘制、频域体绘制和基于蒙特卡罗积分的方法也进行了研究。【19-23】我们注意到，之前的许多工作都是针对患者设置的刚性2D/3D配准进行的，这强调了解剖姿势的精确几何再现，并得益于骨骼特征的高对比度。此外，许多作者使用图形处理单元(GPU)来加速DRR计算。【23】我们的算法使用标准的带有GPU加速的光线投射方法。

【24】DRR算法的基本光线追踪部分只模拟主光子的衰减，而投影图像的形成受到附加物理效应的影响。其中一种因素是散射，这种效应在CBCT投影中有很大的影响。【25】估计CBCT散射的方法有多种，包括波束阻挡阵列、散射核方法和蒙特卡罗模拟。【14】Ruhrnschopf和Klingenbeck对散射估计和校正技术进行了全面的总结。【26】我们的方法是使用模拟CBCT系统的内部蒙特卡罗算法来获得散布估计。与束阻阵列方法相比，这种方法不需要任何特殊设备，更重要的是避免了任何额外的扫描，这会给患者带来额外的剂量。

【27、28】影响投影成像的第二个因素是射束硬化。同样，有很多论文描述了射束硬化和射束硬化校正，但几乎所有的论文都应用于改进CT重建。有两种流行的方法假设受试者完全由水或骨头组成。【29】Yanet等人描述了一种使用迭代方法通过将CT中的每个体素表示为基础材料的和来校正单能量和双能量CT的方法。在他们的DRR算法中，Mooreet等人显式地模拟了多能光束，而不是应用校正。对于我们的DRR算法，我们对后处理中实现的射束硬化校正的分析和测量方法进行了评估。

【30,32】最后一种效果影响是掩盖眩光，即探测器像素向其相邻像素泄漏信号。在我们的机构中，我们测量了在本研究中使用的平板探测器的遮蔽眩光的影响。

我们通过使用高效的射线追踪程序结合散射、光束硬化和遮蔽眩光校正来实现最大化。在本文中，我们提供了一步一步的测量和计算所需的算法。尽管对这些影响的许多描述已经分别进行了描述，但很少有工作将它们合并到一个算法中，以计算真实的drr。

2、 METHOD AND MATERIALS

2.1、定义投影

我们首先定义一个通过物体的CBCT投影
在这里插入图片描述

X_i:探测器像素i的位置(所有坐标都在世界坐标系中)
P_cbct:探测器像素i记录的信号
A₁：管输出系数是否与mAs相对应
$\Phi _{CBCT}\left ( X_{i} \right )$ :以探测器像素为目标的源通量分布
$\alpha_{CBCT}\left ( X_{i} \right )$ :是透射率，也就是光束到达探测器像素时残留的初级粒子的比例
$R_{CBCT}\left ( X_{i} \right )$ :探测器响应函数是否与像素相关
$S_{CBCT}\left ( X_{i} \right )$ :为散射光子对探测器信号的影响

更进一步定义为：
在这里插入图片描述

缺

最后，定义Eq.(1)的一个特殊情况称为泛洪投影是有用的。在光束中没有物体时进行洪水投影，因此在探测器表面上的传输为1，散射为0。这导致了以下简化的数学描述:
在这里插入图片描述

$F_{CBCT}\left ( X_{i} \right )$ :检测器像素i记录的洪水投影信号
洪水投影经验性地捕捉到了源通量分布和探测器响应函数。

2.2、DRR 挑战

对于入射X射线的特定能量光谱，CT图像表示成像对象中材料的X射线衰减系数，相对于水的衰减系数。我们算法的目标是使用一个物体的CT图像来模拟如Eq.(1)所定义的物体的x射线照片(投影)。我们可以将产生投射的物理过程分为主要效应和次要效应。

主要：
成像的主要过程是物体对主光子的衰减。这个过程由对象的LAC分布μ(x)控制，方式如下:
在这里插入图片描述

I₀:光子束的初始强度
I_i：是光束穿过物体后剩余的主光子的强度
从x0到xi进行积分计算

用eq（5）来计算eq（1）中的 $\alpha (_{X_{i}})$ ;
式(5)中的线积分可以用射线追踪程序来估计;
然而，我们进行投影的CT具有无量纲的CT数(相对于水的衰减)单位，而不是绝对的LAC。因此，我们的第一个要求是将CT数转换为绝对LAC，以便以物理上有意义的单位呈现CT。
因此，我们的第一个要求是将CT数转换为绝对LAC，以便以物理上有意义的单位呈现CT。
由于绝对LAC依赖于入射x射线能谱，这一事实使情况变得复杂。因此，必须使用希望模拟其投影的CBCT系统的精确x射线谱来测量通过已知成分物体的x射线衰减。

第二：
我们需要建模的第二个过程是散射、光束硬化和遮蔽眩光。
散射是指被物体偏转但仍被探测器吸收的粒子。
光束硬化是指当多能光束通过衰减介质时，其能量分布向更高能量转移的效应。
遮蔽眩光表示信号从产生它的检测器像素泄漏到它的邻居。这些影响必须从两个方面来解释:

它们必须从用于确定CT数和LAC转换功能的测量中去除
它们必须被折叠到射线追踪程序生成的原始DRR中。

这种分离确保了散射、光束硬化和遮蔽眩光在每个DRR中明确而精确地得到解释，而不是在CT数转换函数中纠缠。

2.3、二级物理效应

2.3.1、光束硬化 Beam hardening

当多能量光束通过材料时，能量较低的光子优先衰减，导致平均光束能量增加。这种效应称为光束硬化。LAC是一个能量依赖的量，因此对于具有统一CT数的给定材料，光束硬化会导致光束在穿过材料时“看到”一个逐渐降低的有效LAC。

为了解释这种效应，我们需要一个函数h(M_bh1)，将实测衰减数据中的硬化累积LAC转换为未硬化累积LAC，以消除用于计算CT数到LAC转换功能的数据中的束硬化。此外，我们需要一个反方向的函数(Mbh0)，以便将光束硬化效应添加到射线追踪程序集成的未硬化累积LAC中。下标“bh1”和“bh0”分别表示有无梁硬化。
我们可以使用相同的数据来实现这两个目标——
一组硬化的累积LAC值，Mbh1，
对应的未硬化的累计LAC值Mbh0，通过不同厚度的测试材料计算。
以mbh0为自变量对数据集进行拟合得到光束硬化函数，以mbh1为自变量对数据集进行拟合得到光束未硬化函数。在本文中，我们提出了基于测量和计算的方法来获得这个数据集。
测量方法的原始数据是一组32个CBCT投影，通过厚度从2毫米到180毫米不等的固体水块。固体水板的位置，使沙发不干扰测量，并有一个最小的体积的板暴露在光束中，以减少分散。我们首先在探测器上选择固体水板阴影下的一个位置作为我们的测量点，并计算每个板厚度的透射率。假设散度可以忽略，可用公式(1)和(4)计算传播，公式如下:
在这里插入图片描述

一般情况下，a0 / a1可以估计为洪水投影mAs与目标投影mAs的比值,或者计算为洪水投影与物体投影的比率在物体投影没有衰减的区域，例如 ,αCBCT(xj)=1.

通过式(2)，可以将每一层固水的透射率转换为硬化的累积LAC，
在这里插入图片描述
在对每个厚度进行这些测量后，我们估计相应的未硬化的累积LAC值如下:

μ0为材料在未硬化光谱中的平均LAC，
l为固体水板的厚度。

虽然我们知道，对于每个平板，我们必须计算μ0。为了做到这一点，我们首先拟合一个8阶多项式函数到我们之前确定的硬化的Mbh1(l)数据，绘制成l的函数，生成一个函数y(l)，如下：
在这里插入图片描述
当厚度l接近0时，我们知道Mbh1(l)和Mbh0(l)将收敛为：

因此，可由式(9)估计多项式函数的线性系数为μ0：

由此，我们可以使用式(8)计算每个板厚的Mbh0(xi)。

我们的第二种方法是通过分析计算的衰减过程，通过测试材料的板，以获得相同的数据集。
这需要三个输入项:

我们的CBCT源的光束频谱，W0(E)
蝴蝶结滤光材料(铝)的LAC分布，μAl(E)
以及测试材料的LAC分布，μmat(E)

通过这些我们计算光束的通过。
我们使用液态水作为试验材料，尽管任何已知衰减系数的材料都是可能的。我们首先让光束通过蝴蝶结过滤器的中间，那里是它最薄的地方，
在这里插入图片描述
出于所有实际目的，这是我们真正的未硬化光谱，因为没有领结滤波器就永远无法获取图像，这是用于校准CT数到LAC转换函数的波束区域，见第二章d。（2.4）
我们现在可以计算我们的测试材料的μ0值，使用
在这里插入图片描述
然后，我们在一个像素一个像素的基础上继续，通过每个像素“看到”的领结滤镜厚度传播光谱。对于探测器像素化，新的光谱计算为：

接下来，每个像素的光谱通过测试材料从0毫米到150毫米的厚度以5毫米为增量进行过滤。对于给定厚度的测试材料j，新过滤的光谱为：
在这里插入图片描述
由此计算出探测器像素i的材料厚度lmat,j对应的累积硬化LAC为：

对应的未硬化累积LAC计算为：

实际应用中，未硬化的累积LAC值无需计算一次，因为它们与检测器像素无关。通过测量或计算获得了我们的数据集Mbh0和相应的Mbh1值，我们准备计算我们的光线硬化和未硬化函数。光线硬化函数(Mbh0)是通过将8阶多项式拟合到绘图的mbh1值作为Mbh0的函数得到的。该功能的输入为未硬化的LAC累计值，输出为加硬化的LAC累计值，具体如下:

在这里插入图片描述

就是未硬化的八阶多项式得到硬化

同样，梁未硬化函数(Mbh1)是通过拟合一个8阶多项式到绘制成Mbh1函数的mbh0值来确定的。该功能的输入为已加固的LAC累计值，输出为未加固的LAC累计值。
在这里插入图片描述
如果用测量方法获得原始的梁硬化数据，每种类型的函数只会产生一个。如果数据是通过计算获得的，每个像素将产生两个函数的唯一版本。

2.3.2 散射 Scatter

散射指的是那些到达探测器的光子之前已经被物体偏转。对于我们的目的，我们必须能够删除和添加分散的效果。为了做到这一点，我们使用了内部的蒙特卡罗(MC)算法，能够模拟只散射投影的对象。蒙特卡罗散射投影可以表示为：
在这里插入图片描述
蒙特卡罗算法也可以模拟洪水投影，

将eq20和eq21相除：

利用我们的CBCT洪水投影，我们将Eq.(22)乘以Eq.(4)得到：

最后，将Eq.(23)乘以Eq.(6)输出因子的倒数，A1/A0,
并且假设 $\Phi _{MC}(x_{i})$ 等价于 $\Phi _{CBCT}(x_{i})$
我们得到了Eq.(3)中描述的散点的估计值，
在这里插入图片描述
这个散射投影可以从另一个投影中添加或减去，以添加或删除散射效果。

2.3.3 遮蔽眩光 Veiling glare

遮蔽眩光是指信号从其原点的探测器像素泄漏到周围的探测器像素。我们使用点扩展函数(PSF)对这种效应进行建模，该函数描述了信号在距离信号源10厘米的区域内的传播。
通过拟合半遮挡投影图像的边缘响应函数，对遮挡眩光PSF进行建模，半遮挡投影图像是通过放置在探测器表面一半以上的铅板获得的。参数函数是高斯和指数项的加权和，分别用于描述遮蔽眩光的短期和长期效应。遮蔽眩光效应可以通过与PSF卷积引入投影，
在这里插入图片描述

下标“vg1”表示有遮蔽眩光存在
“vg0”表示没有遮挡眩光
⊗是卷积运算符。

同样地，通过反卷积PSF可以消除眩光，
在这里插入图片描述

2.4 CT到LAC转换函数

为了计算Eq.(5)中的积分，我们必须将CT数转换为线性衰减系数的单位。这种转换取决于能谱(因此取决于kVp)，在这个能谱上，人们希望模拟CBCT投影。【14】我们的方法与Mooreet等人使用的一般程序相同。我们开始收集由不同材料组成的圆柱形塞。插头是已知尺寸的插入物，取自Gammex电子密度校准模体。
我们首先获取体模的FBCT，以确定每个插头的CT编号。接下来，我们分别捕捉每个插头在给定技术下的CBCT投影。
从这些投影中去除射束硬化、散射和遮蔽眩光后，可以使用公式(2)和(8)计算出每种插头材料的LAC值。
完成此操作后，对作为CT数函数绘制的LAC的二阶多项式拟合将产生我们所需的CT数到LAC的转换函数， $f_{LAC}[CT(x)]$ 其中CT（x）是位置x处的CT编号。

首先，从插头突出部分移除veiling glare
在这里插入图片描述

之后，消除散射点被描述为：
在这里插入图片描述
然后将投影还原为累积LAC：

并应用修正来消除光线硬化

最后，根据累积的未硬化LAC计算LAC

t：plus的厚度7cm
x_c：探测器的中心坐标，在成像过程中，plus的位置

u_{CBCT，plug，vg0sc0bh0}作为CT数的函数，用多项式拟合得到我们的转换函数 $f_{LAC}[CT(x)]$ (如前所述，这种转换功能是特定于x射线管技术的CBCT投影)。为了说明散射、光束硬化和遮蔽眩光校正对最终转换函数的影响，对于未校正和完全校正的数据，用于构建函数的数据点及其多项式拟合如图1所示。

在这里插入图片描述

图一：不同plug材料衰减测量生成的LAC与CT数字数据点和二阶多项式拟合图。衰减测量要么未校正，要么完全校正散射、光束硬化和遮蔽眩光效应。

2.5、DRR工作流

用于产生完全校正的DRR的工作流程类似于用于处理插头凸起的工作流程。

首先，使用x射线体积绘制例程产生原始DRR，我们使用基于光线跟踪的方法来近似公式(5)中的线积分，如下所示：

在这里插入图片描述

K：K是介于x0和xi之间的步数
$\Delta X$ :光线跳跃向量
$\Delta x$ ：射线步长
使用三线性插值在位置X0+k $\Delta X$ 处对CT进行采样。

该算法在GPU平台上的光线跟踪和插值的硬件实现中得到了充分的支持。

在产生原始DRR之后，引入束硬化效应：

洪水投影允许我们再现DRR中CBCT源的通量分布，输出因子允许我们恢复目标投影的束流输出。
然后对该图像进行散射校正：

其中散射投影是使用第二节C.2的方法计算的。
最后，引入面纱眩光效应，对DRR进行了完全校正，

2.6 投影采集

为了验证，通过对象FBCT计算DRR，并与同一对象的实际CBCT投影进行比较。在该研究中测试了两个不同的物体，Catphan^TMCT体模（phantom Laboratory，Salem，NY）和实际患者。在每种情况下，分别在20ms、20 mA和125kVp获得投影。Catphan^TM投影是在全扇模式下采集的，而患者投影是在半扇模式下采集的，每个投影都使用各自的铝制蝴蝶结过滤器。所有图像均在OBI TM CBCT系统(加利福尼亚州帕洛阿尔托的Varian Medical Systems)上采集，该系统使用安装在CLINAC上的G242千伏X射线源和PaxScan 4030CB平板探测器。此外，探测器系统顶部还设置了一个10:1的防散射格栅，由0.5毫米铝盖保护。

3、实验结果

我们通过比较我们的DRR和它们打算匹配的实际CBCT投影来测试DRR算法的准确性。我们计算了两组不同的DRR-一组仅使用测量的通过介质的线性衰减(未校正)，第二组使用Full算法，对LAC测量和DRR投影中的光束硬化、散射和遮蔽眩光进行校正。
这显示了更正的相对重要性。我们使用两个不同的对象作为测试用例-Catphan TMCT幻影和实际患者。
对于这两组DRR，在对原始DRR进行后处理时，使用CBCT洪水投影来测量CBCT系统的输出和源通量分布。
我们强调，这些DRR是检测到的信号强度的绝对计算，没有任何归一化因子或其他数值调整来匹配CBCT预测。
对于定性验证，我们专注于视觉比较，像素强度分布的比较和图像强度直方图的比较。为了进行定量验证，我们将DRR像素强度与实际投影像素强度相关联。

3.1 可视化比较

图2显示了通过CT模体投射的投影，图3和图4显示了从不同角度投射通过实际患者的投射。每幅图中的图像被设置为相同的窗宽窗位以供比较。在图2中，每个DRR的底部都被Catphan ^TMFBCT中的金属伪影污染。这对我们方法的有效性没有影响，因此被忽略了。
在这里插入图片描述

图2：通过Catphan^TMCT体模以0的投影角投影◦. (a）未校正的DRR，（b）完全校正的DRR，（c）实际的CBCT投影。（c）中的红线表示线轮廓采样的位置

在这里插入图片描述

图3：投影角度为92.67◦。(A)未校正的DRR，(B)完全校正的DRR，以及©实际的CBCT投影。©中的红线表示对线条轮廓进行采样的位置

在这里插入图片描述

图4：通过病人的投影角度为178.93◦。(A)未校正的DRR，(B)完全校正的DRR，以及©实际的CBCT投影。©中的红线表示对线条轮廓进行采样的位置。

该评估的目的不是像放射科医生可能提供的那样，提供DRR图像质量的专家分析，而是对我们的算法的输出进行一些基本的视觉观察。我们首先注意到，未校正和校正后的DRR都与CBCT的投影非常相似，最明显的特征都清楚地再现了。通常，未校正和校正后的DRR看起来相似；由于FBCT相对于探测器像素大小的粗略分辨率，两者看起来都不如CBCT投影清晰。这种影响会导致校正的和未校正的DRR中的一些细节丢失，例如图3中的脊柱和椎间盘。

然而，尤其是在患者的投影中，相对于实际的CBCT投影，校正后的DRR显示出比未校正的DRR更真实的对比度。在图3(A)(未校正的DRR)中，患者右上方基本上是饱和的黑色，而实际投影(图3©)显示出相当大的对比度。在校正后的DRR图3(B)中，一些丢失的对比度已经明显恢复。类似地，图4示出了区域，特别是在肺脊上方的图像顶部，其中未校正的DRR饱和为黑色，而校正的DRR和实际的CBCT投影两者都保持可见的对比度。