Fast and flexible X-ray tomography using the ASTRA toolbox

使用ASTRA工具箱进行快速灵活的X射线断层扫描

论文链接：http://dX.doi.org/10.1364/OE.24.025129

项目链接：https://astra-toolboX.com/indeX.html

Abstract

从一系列投影图像中重建物体，如在计算机断层扫描(CT)中，是许多不同应用领域的流行工具。现有的商业软件通常为最终用户提供足够精确和方便使用的重建工具。然而，在使用非标准采集协议或需要高级重建方法的应用中，标准软件工具通常无法计算准确的重建图像。本文介绍ASTRA工具箱。ASTRA工具箱针对多个层析成像应用领域的研究人员，为层析成像投影和重建提供了一套高效、高度灵活的开源工具。讨论了ASTRA工具箱的主要功能，并给出了几个用例。

1. Introduction

X射线计算机断层扫描(CT)的使用并不局限于它在医学诊断中最广为人知的应用。实际上，它查看某些对象内部的能力也使许多其他应用领域受益。例子包括生物医学研究，通常依靠微型ct扫描装置来研究某些物质或药物对小动物的影响[1]；在工业应用中，层析成像重建可以成为conveyor belt上(半)自动化过程质量控制的宝贵工具[2,3]；化学，可用于研究和优化反应及其结果(例如， polyurethane foams)；以及各种成像方式的材料科学，如微型CT、同步辐射设备和电子显微镜[4]。虽然所有这些应用程序共享相同的基本概念(从不同方向测量多个投影图像)，但扫描系统的具体设置可能因应用程序而异。

NRecon用于布鲁克微型ct系统，PhoeniX datosjX用于GE系统，Inspect3D用于FEI电子显微镜)。也有一些商业产品没有与特定的扫描系统紧密耦合，但价格较高（例如，VG Studio MaX，Octopus Reconstruction，Ikeda’s TomoShop）。这些商业软件包的主要优点是由于其庞大的用户基础和用户友好的图形用户界面而具有可靠性。因此，对于大多数CT用户来说，这些工具已经足够了。然而，这些包提供的实现只提供了有限的灵活性，限制了一些用例：

• 重建算法的选择通常限于滤波反向投影(FBP)、Feldkamp-Davis-Kress (FDK)，在某些情况下还包括同步迭代重建技术(SIRT)[5]。然而，其他重建方法(如FISTA)可能会提供更好的结果，特别是在低统计量(高噪声)投影数据的情况下[6]。此外，扫描对象的先验信息不容易整合，无法提高重建质量。这种先验知识在处理有限数据(例如，低数量的投影，有限的角度范围，截断的投影数据等)时特别有用[7]。
• 在大多数CT系统中，X射线源和探测器在垂直于旋转轴的平面上围绕物体遵循预定义的圆形轨迹。然而，在许多应用中，X射线源和探测器阵列必须遵循非常规的轨迹，这意味着它们的投影数据不容易使用传统软件重建。例如，螺旋轨迹现在不仅用于医疗系统，也用于微ct系统的岩石分析[8]。此外，在典型的断层合成装置中，探测器面板位于患者下方，X射线源遵循线性轨迹。同样，在层析中，源和检测器分别固定在薄物体的上方和下方，然后对薄物体进行倾斜和旋转[9]。因此，轨迹是圆形的，但不垂直于旋转轴。此外，在工业过程的质量控制中，只能从几个任意角度获得投影，即使可以在圆形轨迹上从所有方向获取投影，在扫描一个方向上拉长的物体时，改变放大倍率也是有益的[10]。
• 实现的算法的计算效率可能有限，导致重建时间长，特别是对于大型数据集。随着重建框架的计算能力不断提高，X射线探测器面板的尺寸和分辨率也在不断提高。

对于非常规CT系统的用户来说，商业软件包的限制可能会严重阻碍新概念和想法的原型设计。

表1简要概述了最近更新的几种用于层析重建的开源软件包：pyHST2[11]、etattention[12]和IRT[13]。虽然这些包中的每一个都提供了最常见的层析例程的优秀实现，但可以注意到，每个包都专注于特定的目标应用程序。据作者所知，目前还没有软件包比它更适合开发可用于多个应用领域、使用非常规几何和各种约束的高级算法。

本文介绍的ASTRA工具箱(All Scales tomography Reconstruction Antwerp)是一个软件平台，用于解决对灵活、高效和易于使用的层析成像算法开发平台的需求。自2010年以来，它一直由比利时安特卫普大学视觉实验室开发，并于2014年由视觉实验室和荷兰阿姆斯特丹的Centrum Wiskunde & Informatica (CWI)联合开发。它在GPLv3许可下作为开源软件免费提供[14]。为了达到最佳效率，它的大部分代码都是用C++编写的，其核心计算卸载到使用CUDA语言的GPU卡上。无法访问GPU工作站的用户可以求助于CPU代码，使用OpenMP进行加速。该工具箱可以通过MATLAB和Python接口访问，为算法和应用程序原型提供了一个快速而强大的平台。

ASTRA工具箱主要用于研究(例如，应用程序原型和算法开发)，较少用于最终用户的日常使用。然而，由于其开放性，ASTRA工具箱可以与现有的实用软件包相结合。其中一个例子是将其纳入TomoPy[15]，将TomoPy的用户友好性和高效的前后处理方法与ASTRA工具箱的灵活性和优化的重建方法结合在一起，形成一个单一的框架。此外，对于电子断层扫描领域的应用，ASTRA工具箱可以通过FEI的Inspect3D软件获得。

ASTRA工具箱的主要优势在于其固有的灵活性，它有三种形式：

• 凭借其广泛的支持，支持各种非常规采集几何形状。
• 凭借其模块化设计，工具箱的每个小构建块(例如，投影操作)可以独立使用，以便它们可以组合在其他现有的代码库中。此外，ASTRA工具箱可以用于多个应用领域，并且使用基本ASTRA构建块为特定应用程序制定的高级重建方法的实现也可以轻松地应用于其他应用领域。
• 具有可扩展性，可以在局部计算平台上快速评估小型问题的重构方法及其参数。然后，生成的重建策略可以相对较少地部署在大型高分辨率数据集的大规模计算集群上。

本文概述了ASTRA工具箱、它的设计和一些用例。首先，在第(2)节中，介绍了工具箱的关键组件。然后，第(3)节演示了它们如何在传统的锥梁设置中使用。接下来的两节将深入研究允许灵活使用该工具箱的一些关键高级特性：第(4)节关注非常规扫描设置的规范，第(5)节讨论如何轻松地将ASTRA工具箱集成到外部代码库中。最后，第(6)节对本文进行了总结。

2. Framework concepts

ASTRA工具箱中有七个不同的概念和模块(图1)。体积和投影数据对象(第(2.1)节)链接到它们相应的几何信息(第(2.2)节)，两者在投影算法构建块(第(2.3)节)中链接在一起，随后与重建算法对象(第(2.4)节)中的数据一起使用。

2.1. 投影和体积数据

ASTRA工具箱的数据模块存储和管理所有数据集，包括体积数据和投影数据。为了对这些数据执行操作(例如，投影或重建)，工具箱首先需要访问它，如果数据仅存储在MATLAB或Python环境中，则不是这种情况。因此，首先必须分配内存，并且必须将数据从接口层复制到实际的工具箱中。然后，该操作返回一个惟一标识符，该标识符随后可用于引用该特定数据集，这与许多编程语言中的文件输入/输出工作方式非常相似。

2.2. 几何

通过将数据对象链接到相应的体积几何或投影几何来指定扫描系统的实验设置。体几何描述重构对象在其上表示的像素或体素网格。体素通常是各向同性的，但也可以指定为各向异性。这在重建的有效分辨率在每个维度上不相同的情况下是有用的(例如，断层合成和层析)。

投影几何定义了X射线源和探测器相对于体积几何的轨迹。这涉及到探测器像素的数量和大小，以及X射线源和探测器绕物体旋转时的轨迹规格。可用的标准投影几何图形是用于二维数据的平行束和扇形束，以及用于三维数据的平行束和锥形束。

在ASTRA工具箱中指定投影几何可以通过两种方式完成。在其默认形式中，源和检测器遵循围绕体积坐标系原点的圆形轨迹，检测器平面垂直于通过源和原点的直线(并完全对齐)。因此，对于圆锥形梁几何，有五组参数需要指定：检测器像素数量、检测器像素大小、从源到旋转中心的距离、从旋转中心到检测器中心的距离以及投影方向的列表（图2(a)）。使用这种方法，可以指定最广泛使用的投影设置，但非常规的非圆轨迹不能。此外，如果不借助投影域中的插值，就不可能对常见的小探测器不对准进行建模，例如探测器倾斜或探测器移位(当探测器的中心不与通过源和旋转中心的线对齐时)，这一过程将在重建图像中引入不准确性，因此应始终避免。

而不是指定所有这些参数一般固定的对象源距离和探测器大小，而不是假设x射线源很好地围绕旋转中心旋转，ASTRA工具箱还允许用户分别指定每个投影。这被称为矢量几何(图2(b))。以锥束投影模型为例，每个投影图像由体坐标系中的四个三维向量定义：

• 源的位置，$S=(s_x;s_y;s_z)$。
• 探测器面板中心，$D=(d_x;d_y;d_z)$。
• 探测平面的主轴，$U=(u_x;u_y;u_z)$和$V=(v_x;v_y;v_z)$。这些向量的长度对应于每个检测器像素的大小。

从用户的角度来看，指定像这样的投影几何形状更加困难，因为必须仔细指定源和检测器轨迹的计算。但是，它确实允许更高级的设置。此外，探测器对准可以在重建过程中进行建模，而不需要在投影域中进行插值。例如，探测器倾斜可以通过在每个投影中围绕线SD相应地旋转U和V来建模;探测器位移可以通过平移探测器平面(由U和V指定)中的D来建模。

2.3. 投影操作

任何重建技术中最重要的组成部分是投影操作。考虑层析成像的代数形式，正演投影算子然后通过计算稀疏矩阵向量乘法(SpMV)来模拟来自2D或3D体积的投影图像：
$$p:=W\nu ,\text{\hspace{1em}(1)}$$
其中$p$是表示2D或3D投影数据的向量，$v$是表示2D或3D体积数据的向量，W是投影矩阵，由体积几何和投影几何共同定义。反向投影操作通过一个重建体对投影数据进行“涂片”，计算如下：
$$\nu :=W^{T}p.\text{\hspace{1em}(2)}$$
ASTRA工具箱中的投影构建块用于计算W的值，并用于计算投影和反向投影图像。由于这些操作通常占用了层析成像算法的大部分计算时间和能力，因此为了使这些实现尽可能地提高计算效率，需要付出很多努力。为了实现这一点，ASTRA工具箱通过NVIDIA CUDA语言使用GPU卡提供的计算能力[16]。如果用户没有可用的GPU卡，还提供了这些构建块的2D版本的OpenMP加速CPU实现。

2.4. 重建算法

所有层析重建算法都建立在基本的投影和反投影模块之上。这些模块是由算法内部进行配置的。ASTRA工具箱中可用的重建方法包括用于平行和扇形光束数据集的分析重建方法FBP和用于锥形光束数据集的FDK，以及几种代数重建方法，如SIRT, SART和Krylov子空间方法CGLS。所有这些都可以用于2D和3D数据集，并且还允许一些可选的自定义，例如在特定掩码中限制像素的重建，或者在迭代期间强制执行最小/最大约束。算法的配置可以在MATLAB或Python接口中通过将其链接到正确的数据标识符和设置一些特定于算法的选项来完成。重要的是要注意，正向和反向投影构建块本身也可以从外部接口中直接访问。这意味着用户可以很容易地开发新的层析算法或原型，其中很大一部分计算负担被卸载到更快的GPU卡上。

3. Conventional cone beam example

本节演示ASTRA工具箱的界面，供希望对一组常规锥束投影图像进行简单的SIRT重建的用户使用。所提供的代码示例利用了ASTRA工具箱的MATLAB接口。Python接口使用相同的设计和接口原则，只是语法不同。

假设投影数据已经在MATLAB环境中以3D数组的形式存储。由于每台CT机使用不同的文件格式来存储投影图像和几何信息，ASTRA工具箱不提供特定文件格式的读取方法，而是依赖于用户能够做到这一点。

1. 首先，定义了投影几何和体积几何。这里考虑600×600×300各向同性体素网格，每个体素的大小为1 × 1 × 1mm。这个体量的中心位于它的中心切片，在旋转轴上。投影几何图形由图2(a)所示参数指定。每个投影记录为512 × 512图像，像素大小为0.7 × 0.7mm。源到目标距离为1000mm，目标到探测器距离为1500mm。总的来说，360投影图像定义在一个完整的角度范围[0,2π)。

   vol_geom = astra_create_vol_geom(600,600,300); 
   proj_geom = astra_create_proj_geom('cone', 0.7,0.7,512,512, ...
                                           linspace(0,2 * pi,360), 1000,1500);
   

2. 投影数据被复制到工具箱内存中，并分配空间存储重建数据。这些操作的结果是在工作流的后续步骤中使用的两个数据标识符。每个数据对象都链接到相应的几何对象。

   vol_id = astra_mex_data3d('create', '-vol',  vol_geom, 0); 
   proj_id = astra_mex_data3d('create', '-proj3d', proj_qeom, p); % ’p' contains theprojection data
   

3. 有了工具箱可访问的数据，就可以配置重构算法对象。这里考虑一个用于3D数据问题的CUDA加速SIRT实现。该算法在内部配置它将使用的投影构建块。对于这种特殊情况，还应用了非负性约束。这是一种简单的先验知识形式，可以提高重建质量。此配置的结果是一个引用算法对象的标识符。

    cfg = astra_struct('SIRT3D_CUDA');
    cfg.ProjectionDataId = proj_id;
    cfg.ReconstructionDataId = vol_id;
    cfg.option.MinConstraint = 0;
    alg_id = astra_mex_algorithm ('create', cfg);
   

4. 算法标识符用于执行100次迭代。

   astra_mex_algorithm('iterate', alg_id, 100);
   

5. 最后，将重构数据检索到MATLAB存储器中，为后续分析做准备。

   reconstruction = astra_mex_data3d('get', vol_id);
   

4. Non-conventional geometries

本节探讨了非常规投影几何CT的几种不同应用。这些都可以在ASTRA工具箱中使用第2.2节和图2(b)中描述的基于向量的方法进行规范。具体来说，以投影图像的总数为1，构造一个1 × 12的矩阵，其中每一行表示单个投影的几何形状。前三列表示x射线源的位置(S)，第四列至第六列表示探测器阵列中心的位置(D)，第七列至第九列和第十列至第十二列分别表示单位探测器向量(U和V)。在构建锥形束投影几何形状时，将该矩阵和探测器数量提供给工具箱(表2)。

4.1. 旋转计算层析术

在传统的锥束CT中，目标在圆形轨道上成像。然而，对于大型平面物体，如电路板和绘画，这是不可能的，因为它们可能不适合在源和探测器之间的纵向物理位置，或者因为该方向的总x射线衰减可能太高。在这些情况下，旋转计算机层析成像可以通过倾斜物体和旋转工作台，使旋转轴不再垂直于点源和探测器中心排所定义的平面(图3)提供解决方案[9]。这相当于一个固定物体的设置，其中源和检测器遵循平行于z = 0平面的圆形轨迹。这些圆的半径取决于物体倾斜的角度α，以及物体到源或物体到探测器的距离。α的选择也代表了小视场(FOV)(小α)和x射线束与物体的长相交之间的权衡，因此可能有噪声的投影统计(大α)。对于α = 90◦，这种设置相当于一个传统的锥形束设置。

在ASTRA工具箱中表示这种类型的投影几何相对简单。定义∆t为每个探测器像素的大小，将∆s表示为从源到坐标系中心的距离，将∆d表示为从探测器阵列中心到坐标系中心的距离。设 i ∈ { 1 , . . . , L } i \in \{1,...,L\} i∈{1,...,L}表示投影指数，${\theta }_i$表示与每个投影相对应的视角。那么完整的投影几何是：
$$\begin{array}{rcl}{S}_i& =& \left(\Delta s\sin \alpha \sin {\theta }_i,-\Delta s\sin \alpha \cos {\theta }_i,-\Delta s\cos \alpha \right),\\ D_i& =& \left(-\Delta d\sin \alpha \sin {\theta }_i,\Delta d\sin \alpha \cos {\theta }_i,-\Delta d\cos \alpha \right),\\ U_i& =& \left(\Delta t\cos {\theta }_i,\Delta t\sin {\theta }_i,0\right),\\ V_i& =& \left(\Delta t\cos \alpha \sin {\theta }_i,-\Delta t\cos \alpha \cos {\theta }_i,\Delta t\sin \alpha \right).\end{array}$$
在图4中，一个例子显示了一个大而平坦的物体的层析重建。图4(a)所示的CAD模型表示一个1000 × 1000mm的电路板，其中有一系列突出边缘的孔。总共有270个800 × 800的投影(图4(b))是使用NVIDIA Optix光线追踪框架对该CAD模型进行了仿真。层叠倾斜设置为α = 50◦，源到目标距离∆s = 750mm，目标到探测器距离∆s = 1500mm。投影图像应用了真实水平的泊松噪声。随后使用ASTRA Toolbox进行300次迭代CGLS重建，其中两个切片如图4©和4(d)所示。很明显，只有一小部分物体进入了视野。

4.2. Tomosynthesis

层析合成是一种主要应用于医学放射学的成像技术，它可以拍摄少量的有限角度投影图像，并且对x射线剂量的最小化至关重要。例子包括乳腺癌筛查，与传统乳房x光检查相比，断层合成具有更高的准确性[17];以及为重症监护病房的患者提供便携式床上放射照相，这些患者无法在CT或MRI扫描仪中定位。断层合成也可以应用于非医疗用途，如机场的行李检查[18]。

在断层合成中，通过在物体上方沿线性或圆形轨迹移动的x射线源和固定的平板探测器获得投影。图5(a)示出一个示例。在这种几何结构中，只能覆盖有限的角度，导致深度分辨率低于平面内分辨率。因此，重构体通常由各向异性体素组成，其体素大小在z方向上比在x和y方向上大几倍。

设 i ∈ { 1 , . . . , L } i \in \{1,...,L\} i∈{1,...,L}表示投影指数，θi表示每个投影对应的采集角度。探测器平面位于物体中心下方的固定位置，距离∆d。完整的投影几何定义为：
$$\begin{array}{cc}{S}_i=& (\Delta s\sin {\theta }_i,0,\Delta s\cos {\theta }_i)\\ D_i=& (0,0,-\Delta d)\\ U_i=& (1,0,0)\\ V_i=& (0,1,0)\end{array}$$
考虑图5(b)中显示的例子。共获取了45幅人手的层析投影图像。投影角θi均匀分布在[−20◦,20◦)的角度范围内。物体中心放置在探测器面板上方∆d = 65mm处，物体到探测器的距离∆s = 1039mm。检测器计数2208 × 2668像素，尺寸为0.160µm。图像在48kV和0.32 ma下获取。在包含552 × 667 × 30各向异性体素的体积上使用100次SIRT迭代重建这些投影图像(图5©)，体素大小为(0.640,0.640,4.0)毫米。

4.3. 传送带断层扫描

在工业应用中，为了实现利润最大化和客户满意度，精确的质量控制往往是至关重要的。使用计算机断层扫描进行无损检测是实现这一目标的一个有价值的工具[19]。在制造中，CT可用于焊缝[20]、3D打印物体[21]、印刷电路板等的质量检测。在食品加工业中，利用CT可以对不符合一定要求的食品(如内部早腐烂的苹果)在运输前进行丢弃[3]，对肉类(如猪胴体)的切割在知道肉、脂肪和骨头的位置后进行优化[2]。

对于大规模生产环境(生产线)中物体的快速检测，使用传统CT并不方便。在CT扫描仪上单独放置每个物体既昂贵又耗时。为了不影响生产线的效率，可以使用联机层析成像。在这里，一个物体在传送带上移动，在传送带的一侧放置一个固定的x射线源，在它的对面放置一个固定的平板探测器。

当物体以速度v在传送带上平移时，穿过源和探测器的视场，可记录多个锥束投影图像。为了允许在整个角度范围内进行投影，还考虑物体绕其中心轴以垂直于输送带的速度旋$\omega$(图6)。在固定坐标系< x,y,z >，这个几何图形定义为：
$$S=\left(0,-\Delta s,0\right),D=\left(0,\Delta d,0\right),U=\left(\Delta u,0,0\right),V=\left(0,0,\Delta u\right).$$
而不是平移和旋转对象，同时保持源/检测器系统固定，同样的设置也可以从一个固定的对象的角度来描述。在这种情况下，源和检测器围绕对象平移和旋转。固定系统的原点位于物体路径和连接源和探测器中心的线的交叉处的中心。对于每一个时间点i，即物体在输送带上的每一个位置，设$t_i$表示自第一次投影以来经过的时间，定义$R_i$为旋转矩阵
$$R_i=\left(\begin{array}{ccccc}\cos (\omega t_i)& & -\sin (\omega t_i)& & 0\\ \sin (\omega t_i)& & \cos (\omega t_i)& & 0\\ 0& & 0& & 1\end{array}\right).$$
虚投影几何在固定物体坐标系下的函数$<x^{\prime },y^{\prime },z^{\prime }>$则可以表示为：
$$\begin{array}{rcl}{S}_i^{\prime }& =& R_{i}S-(A_0+\nu t_i)R_{i}U,\\ D_i^{\prime }& =& R_{i}D-(A_0+\nu t_i)R_{i}U,\\ U_i^{\prime }& =& R_{i}U,\\ V_i^{\prime }& =& R_{i}V,\end{array}$$
式中$A_0$是在第一次投影时< x,y,z >和< x’ ,y’,z’>坐标系原点之间的距离。

在图7中，显示了一个玩具的内联重建的示例。内联投影数据来源于常规圆形CT扫描的493个投影(图7(a))。换算后，每隔10ms得到100 524 × 3000个投影，源物距离∆s = 121mm，物体到探测器距离∆d = 40mm，探测器像素尺寸∆u = 35.36µm，传送带速度v = 100µm/ms，物体旋转速度为$\omega$=−0.36 =ms, $A_0$ = -50mm。图7(b)显示了在没有物体在传送带上旋转(即$\omega$= 0^◦ms)的情况下中心切片的重建。图7©显示了相同的重建，但有旋转$\omega$=−0.36◦ms。很明显，需要额外的旋转来获得准确的重建。

4.4. 自适应缩放

在锥形束几何中，源到探测器的距离有时被称为缩放因子，因为源离物体越近，它在探测器上显示的越大。理想情况下，选择这样一个物体的投影正好填满探测器阵列。这样，投影图像具有最高的可实现分辨率，而不会引起对象截断和感兴趣区域伪影。然而，考虑一个拉长的物体，如图8所示。在传统的设置中，由于物理限制，源到目标的距离应该根据物体的最长轴与源到探测器的线对齐的方向来选择。这保证了该方向上的最佳分辨率(图8(a))，但导致检测器在其他投影上的次优使用(图8(b))。

使用自适应缩放方法，其想法是在扫描时将物体的位置移动到离源更近或更远的地方。最优的方法是通过扫描物体的凸包预先计算轨迹[10]。这种方法可以获得更多的物体信息，从而提高重建质量。

为了模拟具有可变源和检测器位置的层析成像系统，在桌面微型CT系统SkyScan-1172 (Bruker-MicroCT，比利时)上进行了实验。用直径为7毫米，长度为15毫米的铅笔作为拉长的物体。获得7个全角度数据集，每个数据集包含600个投影，大小为880 × 666，源物距离范围为80.77 ~ 117.01mm。源到探测器的距离固定为216.392mm。

图9(a)显示了最大源物距离下的700次迭代SIRT重构。在此基础上，得到了物体的凸包。对于每个投影角度，根据[10]中描述的技术计算最接近的可能的源位置，并从大于或等于计算的源位置到旋转中心的距离的最小距离获得的数据集中选择一个投影。得到的轨迹如图9©所示，在图9©中可以看到，对于物体具有窄阴影的角度，缩放因子增加了。相应的重建如图9(b)所示，其中自适应缩放方法使对比度略有提高。

4.5. 自动投影对准

以矢量形式考虑投影几何的另一个优点是，当涉及到校正投影不对准时，它具有固有的灵活性。这种错位可能以多种形式出现。例如，旋转中心可能不会正好投影到探测器阵列的中心，并且探测器平面可能不会正好垂直于连接源和旋转中心的线。此外，旋转阶段可能有较小的误差，焦点光斑大小可能不是恒定的，扫描对象在扫描过程中可能会移动。所有这些偏差都可以通过改变投影几何的向量来简单地建模。然而，要做到这一点，首先需要确切地知道在扫描过程中发生了哪些错位，以及错位的程度。最近，有研究表明，通过优化投影几何中的向量，可以自动找到这些不对准参数，从而使重建图像的测量数据与模拟数据之间的投影差最小[22]。

在这里，我们考虑该工作的简化版本，只关注失调因素的一个子集：探测器偏移δU和δV，以及探测器倾斜φ(图10(a))。这些错位通常会导致双棱，如图11(b)所示。当偏移和倾斜因子已知时，可以首先平移和旋转投影图像，然后继续使用传统几何图形进行重建。然而，与其改变数据来适应几何形状，这需要在投影域中进行插值，不如改变几何形状来适应数据。设$R_{i,\phi }$表示以角φ绕向量$D_i-S_i$旋转的矩阵。校正后的投影几何可以表示为：
$$\begin{array}{rlrl}& {\mathbf{S}}_i^{\prime }& & ={\mathbf{s}}_i,\\ & {\mathbf{D}}_i^{\prime }& & ={\mathbf{D}}_i+\delta U{\mathbf{U}}_i+\delta V{\mathbf{V}}_i,\\ & U_i^{\prime }& & =R_{i,\varphi }{U}_i,\\ & V_i^{\prime }& & =R_{i,\varphi }{V}_i.\end{array}$$
有时，偏移和倾斜因子是由扫描系统测量的，可以在与数据相对应的日志文件中找到。然而，在其他情况下，这些信息是缺乏或不准确的。手工校正是乏味的，因此应该避免。图10(b)给出了一种自动估计未知对准参数的算法示意图。在测量投影的基础上，采用常规投影几何(即φ = δU = δV = 0)计算第一次重建。从该重建体中创建单个正演投影图像，并使用配准程序找到最适合测量和模拟投影数据的平移和旋转参数。由于原始重建是使用不对齐的投影几何图形创建的，因此它可能包含一些不对齐的伪影，这也可能影响到registration过程。

为了演示这种方法，请参考图11。在意大利的里雅斯特Elettra的symep波束线上，在TomoLab微型CT上扫描了一个装有糖果的塑料盒和一个高衰减的QRM分辨率幻影。总共获得了1444个890 × 1336的投影，这里只考虑其中的288个。源物距离为6000mm，物探测器距离为1000mm，探测器尺寸为0.375µm。单个投影图像如图11(a)所示。图11(b)所示的未校正重建清晰地显示出严重的不对准伪影。在上述自动对齐方法的4次迭代之后，这些工件已经消失了(图11©)。图11(e)和图11(f)显示了相应的投影差。

5. Algorithm prototyping

假设用户有一个自定义的MATLAB库，其中实现了一些特殊的高级重构算法。一个例子可能是某种形式的总变异最小化(TVmin)，例如使用Chambolle-Pock[23]，目前在ASTRA框架内没有本机实现。定制库通常是为一般用途而构建的，它们的正向和反向投影操作可能作为稀疏矩阵向量乘法(SpMV)执行。为了使用这些库，用户需要负责投影矩阵的构造。虽然ASTRA Toolbox提供了将任何2D几何设置转换为稀疏投影矩阵的功能(参见表3)，但这只能用于非常小的问题。对于实际的数据大小，投影矩阵太大，无法装入主内存空间。此外，由于层析重建的性能是延迟限制而不是计算限制[24]，因此最好总是动态地重新计算投影矩阵，即在每次正向或反向投影期间。ASTRA工具箱为这些投影运算符提供了非常有效的实现，利用了GPU硬件的强大功能。因此，理想情况下，ASTRA工具箱的投影操作可以在自定义重建库的高级逻辑中使用。

最简单的方法是打开外部库，用对ASTRA工具箱中适当元素的函数调用替换所有SpMV出现的情况(表4)。虽然可能，但这是不可取的，因为用户可能不熟悉自定义库的内部代码，并且维护其代码可能相当困难。

更好的方法是使用最近的MATLAB库Spot-Toolbox[25]。这允许使用自定义逻辑覆盖像矩阵乘法这样的操作。ASTRA Toolbox for MATLAB与这样一个称为opTomo的点算子相结合[26]。该算子可以配置设置的投影和体积几何，之后MATLAB将在许多方面将其视为正常矩阵。因此，它可以作为普通函数参数传递给外部库。但是，每次MATLAB检测到该矩阵与向量的乘法或转置矩阵与向量的乘法时，分别调用ASTRA正投影或反投影代码(表5)。

作为一个例子，考虑同步迭代重建技术(SIRT)。这个代数重构法将重构问题视为对系统的求解：
$$W_{\nu }=p,\text{\hspace{1em}(3)}$$
其中p是包含投影数据的向量，v是表示未知重构体积的向量，W是将体积几何映射到投影几何上的投影矩阵。由于公式(1)无法直接求解，SIRT采用如下迭代方案：
$${\nu }^{(k+1)}={\nu }^{(k)}+C{W}^{T}R(p-W{\nu }^{(k)}),\text{\hspace{1em}(4)}$$
其中C是包含W的列和逆的对角矩阵(即$c_{jj}=1/{\sum }_i{w}_{ij}$)， R是包含W的行和逆的对角矩阵(即$r_{ii}=1/{\sum }_j{w}_{ij}$)。

在图3中，给出了一个MATLAB函数，该函数在给定投影矩阵(W)、投影数据§、初始估计(v⁽⁰⁾)和若干迭代的情况下执行SIRT重构。在实践中，这个脚本并不有用，因为SIRT在ASTRA工具箱中也有一个本机实现，它要快得多。然而，这确实提供了快速创建新想法原型的灵活性。要使用SIRT函数，需要从ASTRA工具箱中提取一个投影矩阵，并将其作为参数传递给该函数。这种方法需要大量的内存，并且没有利用ASTRA工具箱的高效实现。通过简单地链接ASTRA工具箱和表4中的脚本，使用ASTRA工具箱实现在GPU上运行投影操作，但是代码从6行变成了10行，可读性差了很多。还要注意，函数参数列表略有不同，这可能会破坏其他代码段。通过使用Spot工具，如表5所示，将两者的优点结合在一起。

在性能方面，图12(a)显示了这些不同的SIRT实现在128 × 128体积上具有180个投影的100次迭代的时间，在包含运行在2.60GHz的英特尔酷睿i7-6700HQ CPU, NVIDIA GeForce GTX 960M GPU和32GB系统内存的系统上执行。在本节介绍的三种实现中，原始脚本(图3)是最慢的，因为它没有利用有效实现的ASTRA Toolbox投影操作。naive(表4)和opTomo(表5)代码都要快得多，而opTomo代码只是稍微弱一点。图12中还包括SIRT的两个本机ASTRA实现。完全在CPU上运行的那个，比任何这些MATLAB实现都要慢得多。这是由于MATLAB也使用GPU计算来执行SpMV。ASTRA基于GPU的实现比最快的MATLAB代码快了一倍以上。这主要是由于在这些MATLAB版本中不断地来回复制数据。在图12(b)中，使用100次迭代的CGLS重建算法重复相同的实验。

在图13中，100次SIRT迭代的重建次数绘制为体积大小(2D和3D)的函数。在每个模型中，共考虑了180个平行光束投影角度。本实验在两个不同的系统上进行：系统A，如前一段所述；系统B, 12核 Intel Xeon E5-2630计算服务器，运行频率为2.30GHz，配备64GB内存和强大的Tesla K20 GPU卡。

6. Discussion and conclusions

本文介绍了ASTRA工具箱，并介绍了几个用例，展示了它的主要特性。这些包括它的免费，开源的性质[14]，它的能力来描述几乎任何2D和3D投影设置使用矢量几何，以及它的灵活性，关于它可以包含在其他框架的容易。结合起来，它们使ASTRA工具箱非常适合新应用，新几何设置和新重建方法的快速原型。此外，由于其有效实现的构建块，这些原型可以直接升级到实际的数据大小和用法。

然而，应该注意的是，ASTRA工具箱绝不是处理层析成像重建时要考虑的第一和唯一工具。许多用户可能更喜欢商业软件包，因为它们提供易于使用的图形用户界面，这是ASTRA工具箱所不具备的。另外，由于ASTRA工具箱在应用领域、扫描设备和协议方面具有灵活性，它无法提供对各种文件格式的即插即用支持。因此，为了使用工具箱，用户需要了解这些文件格式，以及在MATLAB或Python层中解析和处理它们的技能。此外，ASTRA工具箱中捆绑的算法仅限于重建方法，不包括典型的预处理(例如，平场校正，相位检索)和后处理(例如，分割，形态操作，网格生成)算法。为了实际使用，该工具箱必须与TomoPy等其他工具相关联[27]。这些缺点限制了ASTRA工具箱的目标受众主要是具有计算机科学专业知识的研究人员和用户。

虽然这似乎是一个巨大的问题，但这些研究团体实际上相对较大，而且更多的研究团体出现在文献中。ASTRA工具箱的应用范围很广，包括(i)材料科学，通过电子断层扫描[4]、同步辐射衍射对比断层扫描[28]、微ct系统、层析成像系统[9]等；(ii)生物医学，通过微型CT[1]、断层合成装置等；(iii)工业用途，通过传送带系统[3]、微CT系统等；(iv)安全目的[18]。

综上所述，ASTRA工具箱是一个很好的平台，可以弥合应用科学家与数值数学和线性求解领域研究人员之间的巨大差距[29]。这样，新的和先进的数值求解器可以在实际数据上进行测试，从而使两个社区都受益。