Arbitrarily large iterative tomographic reconstruction on multiple GPUs using the TIGRE toolbox

使用TIGRE工具箱在多个GPU上进行任意大的迭代层析重建

论文链接：https://arxiv.org/abs/1905.03748

项目链接：https://github.com/CERN/TIGRE

Abstract

随着时间的推移，层析图像的大小一直在增加，而计算层析重建的GPU的内存大小也在增加，但它们的速度不够快，无法重建最大的数据集。这个问题通常是通过在大量GPU集群中重建数据来解决的，这些GPU集群中有足够的设备来适应测量的x射线投影和重建的体素。这通常需要数十个GPU，这是一个非常经济昂贵的解决方案。在该领域，使用设计用于仅使用一个或几个GPU进行重建的单节点机器更为常见，但目前的软件不允许迭代重建不适合这些GPU的体素。在这项工作中，我们提出了一种策略，可以在单个节点上使用任意数量的GPU和任意小的内存，有效地执行任意大图像的迭代重建所需的操作。提出了投影和反投影算子的策略，以及两种易于推广到其他投影类型或正则化器的正则化方法。提出的改进还加速了单个或多个GPU上较小图像的重建，为时间紧迫的医疗应用提供了更快的代码。生成的算法已被添加到TIGRE工具箱中，TIGRE工具箱是用于一般CT的迭代重建算法的存储库，但这种节省内存和问题分割策略可以很容易地适用于任何其他基于GPU的CT代码。

1. Introduction

近年来，由于X射线源、探测器技术和计算资源的进步，层析成像应用中的图像尺寸显着增加，同时可用的重建和处理时间要求越来越小。在工业或科学应用的计算机断层扫描(CT)中，可以覆盖2m× 3m探测器区域的扫描仪已经在商业上可用，探测器像素尺寸约为0.2mm，从而产生非常高分辨率的投影，从而产生非常高分辨率的图像。重建这些图像是一个很大的挑战，因为它们的大小可以达到数十到数百GiB，比当前图形处理单元(GPU)的内部存储器还要大，而GPU是CT中最常用的重建设备。另一方面，CT的工业和医疗应用也在推动更快的重建。CT用于制造部件的质量控制和计量学需要快速扫描时间，因为需要增加样品的吞吐量[1] [2]。在医疗应用中，特别是在自适应和图像引导放射治疗中，短的扫描和重建时间是正确靶向肿瘤的基础[3] [4]。

特别是在医疗领域以及其他应用中，迭代算法越来越普遍，因为它可以在低数据质量和数量的场景下提供更高质量的图像，减少扫描次数和辐射剂量[5] [6] [7]。然而，迭代算法明显慢于基于标准滤波的反投影重建方法，不利于快速计算时间和更大图像尺寸的目标。

用并行计算重建大型图像并不是一个新的或未被探索的问题，关于该主题有大量的文献，特别是滤波反投影(FBP)及其锥束修正，即Feldkamp Davis Kress (FDK)算法。对于平行光束情况，通过将问题分割成轴向切片并独立求解[8] [9] [10] [11]，可以很容易地将问题分布在GPU或高性能计算机(HPC)中。

对于锥束层析成像，图像不能被简单地分割成独立的碎片，从而使可能的并行化过程变得复杂。在这种情况下，基于cpu的方法是最常见的并行化方法，因为它们不需要更改代码，只要在一个足够大的基于CPU的计算机上执行，比如HPC集群[12] [13] [14]，它们就可以工作。然而，GPU在CT重建方面具有显著的优势，并且已经提出了几种单GPU和多GPU方法[15] [16] [17]，通常使用消息传递协议(MPI)[18]来允许在每个节点具有几个GPU的大型集群中执行。这种方法可能是大图像CT重建的最佳方法，因为如果优化图像的分区以最小化内存传输，则可以实现高效的计算[19] [20]。

然而，目前可用的文献和软件在多个GPU上重建大型图像时有一个主要缺点：整个投影集和图像必须适合GPU内部RAM。要重建具有2000³个体素的图像体素，这需要大约8个最高内存的GPU，而对于4000³个体素的体素，则需要超过50个GPU，具体取决于投影大小。使用这些资源的计算机并不常见，而且非常昂贵。大多数工业、科学和医疗中心无法访问这种多GPU高性能计算集群，而是使用专用的单节点机器进行CT数据处理，每个机器有几个GPU(最多2到8个GPU)。一些已发布的FDK算法的GPU实现可以通过将图像分割成独立的片段并独立处理每个片段来克服这个问题[21] [22]，但没有在任意小的GPU上进行任意大的迭代重建的通用方法。在这项工作中，我们提出了一种分割投影和反投影操作的策略，以便有效地单CPU多GPU执行大型CT图像。

大多数多GPU层析成像重构代码拆分数学运算符，使标准重构算法保持不变。然而，另一种方法是用数学方法来分割问题[23] [24] [25]，提出不同的分布式优化方法，使用多个局部收敛的独立计算来确保重建的全局收敛。每个线程都可以在单独的机器上执行，节点之间的内存传输有限。这些方法可以在任意大的计算机网络中执行，但使用不同的算法来解决重建问题。

本文的其余部分结构如下。首先，在方法部分，概述了TIGRE和迭代算法，并回顾了GPU术语。本节继续解释投影和反投影操作的建议方法，并以如何加速在GPU上实现的正则化操作的描述结束。然后，我们给出了不同系统的执行时间的结果，这些系统具有不同的图像和检测器尺寸，以及用测量数据重建的两个示例。

2. Methods

高效的多GPU层析重建策略主要基于重新安排计算和内存传输操作，从而使计算量最小化，而内存传输和计算执行之间的重叠最大化。

基于GPU的层析迭代重建工具箱(TIGRE)是一个基于MATLAB和Python的GPU加速层析重建工具箱[26]。它特别集中于提供了十多种迭代算法的实现，以及FBP和FDK算法。

从数学上讲，迭代算法求解线性化模型：
$$Ax=b+\tilde{e}\text{\hspace{1em}(1)}$$
其中$x\in {\mathbb{R}}^{N_e}$是图像元素的向量，$b\in {\mathbb{R}}^{N_p}$是吸收测量的向量，$A$是描述射线和元素之间相交长度的矩阵。文献中有多种使用迭代方法求解公式1的方法，但它们都有一些共同点：每次迭代中至少需要$Ax$(投影)和$A^{T}b$(反投影)操作中的一个。这些是计算量很大的操作，占用了99%以上的执行时间。一些算法还包含其他计算量大的操作，如正则化操作。这些操作是在使用CUDA语言的NVIDIA设备的GPU加速代码的TIGRE中实现的。如果要重建任意大的图像，这些核心操作是需要适应的。分割和收集部分结果的低效策略可能导致不必要的计算和内存传输。本文的目标是针对一般情况优化拆分策略。

TIGRE的架构是模块化的，因此所有的GPU代码都独立于使用它的算法。因此，通过调整GPU代码以能够处理任意大的重建，TIGRE也将自动处理这样的图像。然而，模块化设计有一个缺点。为了通过将计算量大的操作转换为“黑盒”函数来实现设计新算法所需的更高级别抽象，需要在内存管理和算法特定优化方面做出牺牲。具体来说，TIGRE将内存从CPU移动到GPU，反之亦然，在每个GPU调用。这允许新方法的快速原型，但意味着算法可能比完全优化的GPU代码执行得稍微慢一些。该架构还带来了另一个限制，特别是对于多GPU代码：由于内存在高级语言(MATLAB和Python)中被分配和释放，它存储在可分页内存中。简而言之，这是由操作系统(OS)管理的内存，如果操作系统认为有必要，可以拆分、移动或不完全存储在RAM中。另一种选择是页面锁定或固定内存，顾名思义，这些内存被锁定在某个RAM位置。这种类型的内存不仅具有更快的传输速度(因为它可以由GPU管理而无需CPU干预)，而且还允许加速代码的重要功能，异步内存传输。

CPU对GPU指令的调用可以是同步的，也可以是异步的。如果调用是同步的，CPU等待直到指令完成后再继续，而如果是异步的，CPU将继续执行进一步的指令，不管GPU是否终止了它的任务。排队操作，但要确保在满足要求之前没有一个操作被执行(例如，一些计算完成，特定的内存传输结束)是在多GPU系统上最小化总计算时间的关键特性，因为它允许排队同时进行存储器读取或写入以及同时进行存储器管理和计算执行。由于TIGRE依赖Python和MATLAB来分配此内存，因此内存最初将被分页，因此只能同步传输，除非事先显式固定，这是一个代价高昂的操作。虽然前面提到的限制会导致这里提供的代码效率稍低，但也给出了在没有TIGRE限制的情况下如何实现多GPU策略的讨论。

在介绍所使用的加速策略之前，值得一提的是，这项工作关注的是如何组织计算和内存传输，但没有描述(也没有约束)内核(计算函数)本身。文献中有几种算法可用[27] [28] [29] [30]在GPU上有效地实现投影和反投影，代码之间的公平比较是一项具有挑战性的任务，因为代码通常不可用，执行时间是特定于硬件的。我们并没有声称在这项工作中使用的内核是最快的，但是我们提出的多GPU策略适用于大多数(如果不是全部的话)文献中的投影和反投影算法。

2.1. 投影

投影运算$Ax$是用给定的图像估计$x$来模拟理想的X射线衰减。有几种已发表的算法来近似该运算，其中两种算法在TIGRE中实现：插值投影[27]和Siddon方法的射线驱动相交投影[31]。虽然每种方法对X射线路径上的积分计算不同，但内核是使用相同的结构组织和执行的。

内核在GPU上排队执行，使用大小为$N_u\times N_v\times N_{angles}$的3D线程块，其中每个线程处理单个检测器像素的积分，如图2所示。通过在GPU中使用纹理的邻近内存缓存功能，由于更高的内存读取延迟，这种内核结构允许更快的执行。这是因为内存读取需要数百个周期才能完成，但将图像存储在纹理内存中允许3D缓存。因此，读取相同或邻近元素的并发线程会提高缓存命中率，直接转化为更短的内存读取时间和更快的执行速度。除了更快的内存读取之外，纹理内存还允许硬件实现三线性插值。关于投影法的更多信息，读者可参考Palenstijn等人[32]的技术说明。每次投影内核启动时，都有足够的线程块排队，以覆盖投影的整个高度和宽度，但深度只有一个。因此，每次内核启动都会计算$N_{angle}$的整个投影。

然而，为了在锥形束几何形状(如Circular, C-arm或螺旋扫描等)中执行投影操作，需要完整的图像，因为由于X射线的形状，不同的图像分区将对探测器中的重叠区域产生影响。因此，如果图像不适合GPU，则需要为每个图像分区分别计算投影，然后再累积。累积过程非常快，执行时间大约是投影内核启动所需时间的0.01%。

在多个GPU上进行投影操作的一种常见方法是根据图像大小分配投影的内存，然后执行内核。然后进行收集步骤，将重叠的检测器块进行累积。如果图像和投影适合可用的GPU内存，这是一个很好的方法，但对于大型图像，情况并非如此，因此变量需要部分存储在主机RAM中。一旦这成为必要，那么最小化GPU上不必要的内存量是很重要的，以便最大化每个单独映像分区的大小(即最小化分割映像的数量)。

我们提出的方法可以总结如下：只有两个投影内核启动的内存将在每个设备中分配($2\times N_{angles}$单个投影)，其余的GPU内存将用于图像。这两个投影块将充当缓冲区，在内核执行期间将计算的投影复制到CPU时存储它们。当一块用于存储正在计算的投影时，另一块将把内存复制到CPU RAM。这不仅允许我们最小化所需的分割次数，还减少了执行时间，因为不需要额外的时间将结果从GPU复制出来。在多GPU系统中，每个GPU将使用整个图像计算一组独立的投影。如果图像不完全适合每个GPU，图像被分割成相同大小的轴向切片堆栈，尽可能大。此外，将分配一个具有相同大小的额外投影缓冲区。这些将在第一个图像切片被前向投影后使用，因为它们将加载之前计算的部分投影数据集到GPU RAM中，以便在每个投影内核完成时在GPU上积累。

以正确的顺序将这些操作排队是关键，因为内存副本将暂停CPU代码，直到完成，从而阻碍并发操作的执行。一旦图像切片已经在GPU中，首先将异步执行投影的内核启动。然后，如果需要，从先前计算的图像切片的部分投影从CPU复制到每个GPU。复制完成后，CPU异步排队到累积内核，最后将前一次计算的结果复制到CPU。这些指令中的每一个都是在所有可用的GPU上同时执行的，而不是一次在一个GPU上执行所有指令。流程的时间线如图3所示，算法1显示了伪代码。

在一个和两个GPU的情况下，经验测试表明，如果由于缺乏内存而需要对图像进行分区，在进程开始时对图像内存进行页面锁定会导致更快的总体执行速度。由于CPU-GPU传输速度更快(在PCI-e Gen3上从大约4GB/s提高到12GB/s)，即使需要大量时间进行页面锁定也是如此。在两个以上的GPU上，页面锁定总是带来改进，因为它支持同步复制。

2.2. 反投影

反投影操作，$A^{T}b$，沿着X射线路径向源模糊探测器值$b$。TIGRE使用基于体素的反投影，有两个可选的权重，FDK权重和“伪匹配”权重，它们近似于光线驱动的相交投影的伴随线[33]。两个反投影都执行$N_x\times N_y\times N_{angles}$线程块，其中每个线程计算$N_z$体素更新，如文章[34]和[35]中详细描述的，如图4所示。类似于投影运算符，这种执行顺序增加了处理器的占用并优化了缓存命中率，减少了总体计算时间。因此投影需要存储在纹理存储器中以启用3D缓存。

与投影类似，文献中常见的反投影方法要求在计算反投影体素更新之前或期间在GPU RAM上分配整个投影集。这种方法不能最大限度地减少内存使用。

算法2给出了多GPU分割算法。

与投影操作不同，在每个GPU中分配两个大小为$N_{angle}$的缓冲区用于投影。在计算开始之前，图像被分割成大小相等的体素轴向切片堆栈，并在GPU之间分配。如果总图像(加上缓冲区)不适合设备之间的GPU RAM总量，则添加一个图像块队列，每个GPU因此负责向后投影多个图像片。每个GPU使用所有的投影，但在体素更新计算的同时执行CPU-GPU内存传输。除非图像大小远远小于检测器大小(在实际数据集中不会发生这种情况)，否则内存传输应该足够快地完成。执行时间线如图5所示。经验测试表明，当单个GPU需要计算多个图像片段时，对图像内存进行分页锁定比使用可分页内存更快。

2.3. 正则化

层析成像重建算法可以包括用于强制用户定义约束的先验信息。虽然文献中有各种各样的正则化器[36] [37] [38]，但总变分(TV)正则化可以说是CT算法中最常见的[39] [40]。因此，我们在这里关注电视类型约束。然而，用于加速和分割问题的逻辑可以应用于其他形式的基于邻域的正则化器。

TV正则化为图像重建问题增加了一个额外的约束：梯度的稀疏性，强制图像分段平坦。TIGRE中使用了两种不同的最小化TV约束的格式，它们出现在文献中不同的算法中：梯度下降最小化方法和Rudin-Osher-Fatemi (ROF)模型最小化方法[41] [42]。这两种代码都使用不同的最小化算法来减少总变化(TV)，但有两个共同点。首先，两种方法都以图像体素作为输入，返回大小相同的图像体素，即算法不是针对CT的。其次，与其他一些正则化子类似，它们需要相邻体素的知识来更新每个体素值，因此是耦合运算。这使得它们很难在多个设备上并行化，因为每次迭代都需要体素值的同步，需要GPU之间的通信。这两种方法通常需要数十或数百次迭代才能收敛。

虽然这些正则化器的优化是特定于算法的，但拆分方法与图6中描述的相同。这种类型的单体素相邻的一个有趣的特性是，如果在包含相邻图像分割的体素值的图像的边界上创建重叠缓冲区，则缓冲区的深度等于图像块可以执行的独立迭代的数量。这意味着可以运行$N_{in}$数量的正则化函数最小化的独立迭代堆栈，然后用其他GPU上新计算的体素值更新缓冲区值。然而，缓冲区的深度直接转化为更高的计算需求，因为每个GPU都需要在相同的片中重复操作。在我们的工作中，我们发现在每个分割边界中，深度值为$N_{in}=60$在最小化内存传输步骤和计算时间之间具有最佳平衡。此外，算法在每次迭代中可能需要一个范数。如果存在需要整个图像的操作，则需要做出特定于算法的决策。在电视的情况下，每次迭代都需要图像更新的范数。经验测量表明，在每次迭代中不同步并沿图像样本均匀分布逼近范数对算法的收敛性和结果的影响可以忽略不计。

然而，当图像和辅助变量不完全适合跨设备的GPU RAM总量时，这种方法就会受到严重影响。由于正则化优化器通常需要输入图像的多个副本(例如，TIGRE中的ROF最小化器需要5个副本)，因此很容易达到GPU RAM的上限。如果是这种情况，映像分割需要不断地读写到CPU RAM中，这会严重影响整体性能。在这些情况下，为了最大限度地提高数据传输速度，内存被分配并固定在CPU RAM中。在副本适合GPU内存的情况下，缓冲区也被分配并固定在主机中，但大小要小得多。

3. Results

进行了两个实验来评估所提出的多GPU代码。首先，使用不同的问题大小和GPU数量，对投影和反投影操作本身执行详细的实验。分析了在不同情况下的加速和性能。其次，使用多GPU代码和TIGRE工具箱重建不同尺度和不同X射线机的两个大扫描。

3.1. GPU调用性能分析

为了测试多GPU代码的性能，CUDA代码在多达4个GPU上以不同的图像大小执行。两台机器用于实验，规格非常相似。对于1和2 GPU节点，使用一台带有AMD EPYC 7551P处理器(256 GB RAM)和2台NVIDIA GTX 1080Ti (11 GB RAM)的工作站，通过独立的PCI-e Gen3 x16端口连接。对于3和4 GPU测试，使用了Irridis 5 HPC集群上的GPU节点，2个Intel Xeon E5 2680处理器，128 GB RAM和4个NVIDIA GTX 1080 Ti通过独立的PCE-e Gen3 x16端口连接。使用不同的机器仅仅是因为2 GPU工作站上的RAM更大。虽然处理器不同，但经验测量显示，在HPC GPU节点中执行1和2 GPU代码时没有显着的时间差异，因此结果具有可比性。

图7显示了不同大小(N)和GPU数量的总时间(计算加内存传输)。实验使用N³个图像体，每个图像体有N²个检测器像素和N个投影角度。图8显示了与单个GPU代码相比所需时间的百分比。

TIGRE有两种投影和反投影类型，但这里我们只对投影和FDK加权反投影的默认射线-体素相交算法进行计时。替代插值投影不会被任何TIGRE算法调用，因为它比光线体素版本慢。为了完整性，它被包含在工具箱中，并且在测试中给出了几乎相同的结果。具有匹配权值的反投影比FDK加权的反投影慢10%-20%，并且仅在基本需要匹配的反投影时使用，例如使用CGLS或FISTA算法时。除了较慢的内核，多GPU实现是相同的，因此它们的性能结果没有显示。

结果显示了预期的性能增益。计算上，复制GPU的数量应该会使内核时间减半。这在很大程度上是正确的，唯一的区别是由于每个内核中的最后一个计算块的nangle值小于它们对应的值，以及投影积累所需的额外计算时间。没有达到理论速度提高的另一个原因是由于内存管理和传输时间。在较小的问题规模下，内核的速度足够快，以至于时间被GPU属性检查和内存传输所支配。计算内核足够小，可以隐藏多个GPU可以提供的任何改进。对于反投影并行化大小为N = 128的特殊情况，执行速度较慢，但总体时间小于20毫秒。随着N的增大，这两种操作的计算成本都在增加，因此测量到的加速比接近理论值(2、3和4个GPU分别为50%、33%和25%)。这在正向投影中可以看得更清楚。添加更多的GPU通常会使比率更接近理论极限，直到有足够的GPU，每个GPU的计算负载与所需的内存传输相比太少，就像本实验中最小尺寸的情况一样。

由于两个主要原因，反投影的缩放效果不如投影。首先，它更快。这意味着内存管理对总体执行时间的影响要大得多。此外，反投影需要固定图像内存，这增加了显著的开销。在投影操作中，内存在执行之前是保留的，但是操作系统会在写入时分配内存。

作为分割过程的一个例子，对于大小N = 3072，单个GPU节点需要11个图像分区，而2个GPU版本需要6个分区进行反投影。投影只需要10和5个分区。差异是由于投影使用的$N_{angles}$值较小。

由于这项工作提出了GPU加速的内存处理优化，因此测量内存对代码总执行的影响是特别有趣的。图9显示了不同尺寸和GPU数量的投影和反投影的结果。这些操作分为三类。计算包含内核启动的时间，其中包括并发发生时的并发内存拷贝。然后加上固定和解锁CPU内存的时间。在反投影操作中，这个时间占总时间的较大一部分，因为它迫使CPU分配内存，而在投影中，内存已经分配，因为图像体素已经存在。有些问题规模没有这个部分，因为页面锁定是一个缓慢的操作，所以没有执行它，因为它带来的改进少于它的成本。最后，还有其他内存操作。这包括将图像复制进或复制出GPU并且没有发生计算的时间。它还包括内存释放和其他次要的内存操作。

这种分析也有助于更好地解释图9。我们可以看到，在投影中，即使是非常小的图像，计算时间也占据了大部分的执行时间，而在反投影中，即使是512³体素的体素，使用超过1个GPU的计算时间也不到一半。

3.2. CT图像重建

为了展示多GPU加速工具箱的功能，重建了来自两个不同扫描仪和不同分辨率的两幅图像。使用的计算机是前面描述的具有256 GB RAM的双GPU机器。第一次扫描的是一颗烤过的咖啡豆，第二次扫描的是一块鱼龙化石。

Coffee Bean on Zeiss Xradia Versa 510

第一个样品是烤过的咖啡豆。在Zeiss Xradia Versa 510上进行扫描，80kvp，87.5µa (7w)。未使用X射线源过滤。源物距离16mm，几何倍率为9.48，宏观物镜倍率为0.4。使用2000 x 2000像素的检测器测量面板移位扫描，将检测器向左和向右偏移，然后缝合成4000 x 2000投影。每次位移都是在360度以上的6401个等间距旋转角度上获得的，曝光时间为14秒，扫描时间为45小时。利用这些扫描参数，获得了3.653µm的体素尺寸。

选择该扫描图使用CGLS算法进行重建。由于该算法需要对重建图像的大小进行多次复制和投影，因此无法在将所有数据保存在系统RAM中的情况下重建整个4000 × 4000 × 2000体素的数据集。因此，投影被裁剪成900 × 3780大小的投影，只有2134个角度(三分之一)用于重建。这将导致需要29 GB内存的投影集。重建尺寸为12.2 × 12.2 × 3.28 mm³的3340 × 3340 × 900体素图像，需要40GB的内存。如前所述，重构机器上的2个GPU在它们之间有22GB或RAM。

图10显示了30次迭代后使用FDK算法(左)和CGLS算法重建的切片。CGLS重建耗时4小时21分钟。虽然在全角度采样时，FDK算法提供了高质量的图像，但当使用三分之一的角度时，它开始在重建图像上显示噪声伪影。CGLS重建对这种噪声具有更强的鲁棒性。

Ichthyosaur fossil on a Nikon/Metris 225kVp/450kVp custom bay

用于此大型图像重建的数据集是鱼龙化石的一小部分[43]，特别是鳍的尖端。在Nikon/Metris 225 kVp/450 kVp定制托架上，360 kVp，3.37µa进行扫描。源的滤波使用一个4毫米的铜滤波器进行。源与目标的距离为1.564 m，源与探测器的距离为2m。所使用的检测器有2000 x 2000像素，在扫描过程中向样品的左右偏移，在图像融合后得到0.8×0.4 m²的有效检测器。以6400个等间隔角间隔均匀取样。重建图像尺寸为3360 × 900 × 2000，各向同性体素尺寸为0.156 mm。

使用OS-SART算法进行重建，这是经典SART算法的一种变体，使用部分投影子集而不是每个投影更新体素。子集大小为200个投影，使用50次迭代执行，执行时间为6小时40分钟。为了将所有数据保存在系统RAM中，只使用了2000个均匀采样的投影。图像为14.5 GB，投影子集使用了62 GB。图11显示了重建图像的标本和切片。

4. Discussion

结果显示了一个清晰的画面：使用提议的内核启动结构，可以使用多个GPU轻松地并行重建非常大的图像，从而最小化所需的图像分区数量。他们还展示了使用多个GPU如何几乎线性地加速重建，特别是对于较大的图像。但是，可以在本文的结构和代码中添加一些可能的改进。

同步内存传输和计算队列策略是一种改进，它允许多GPU代码在这项工作中以很少的开销实现高并行性，因为在计算之前或之后不需要等待投影被复制到GPU中或从GPU中复制出来。然而，在Zinsser等人[35]的工作中，提出了一种可用于重叠图像存储传输的策略。在他们的工作中，投影完全保存在GPU内存中，但在执行计算时，部分以增量大小复制。类似地，它们在单个GPU中将图像分成两部分，并开始将其中一部分复制到CPU，而另一部分完成计算。每个GPU将图像分割成不同的切片也可以应用于投影操作。这将允许两个操作符稍微减少内存传输时间。但是，为了确保同时进行计算和内存传输，需要将映像拆分为非常特定于计算机和问题的大小。内存传输速度和不同的投影/图像大小是影响图像需要分割多少次的主要变量。由于TIGRE试图尽可能地灵活，这种改进并没有实现，但这里需要指出的是，特定于应用程序的CT算法可能需要考虑它。在处理投影时，不存在内存操作可能比内核花费更长时间的问题，因为内存传输时间总是比内核的计算快得多，除非投影和图像大小之间存在不切实际的巨大差异。

一些改进可以在GPU代码中实现，特别是在反投影中，如果可以给代码一个已经预分配的内存块。由于工具箱的模块化设计，每个内核调用都为其输出分配内存，但是由于迭代算法的性质，很可能已经存在要覆盖内存的变量。我们决定不包括这一改进，并保持TIGRE的模块化设计。

内存的上限由可用的CPU RAM设置，但是GPU代码有其局限性。如果考虑到内核最快执行的值(投影操作中的$N_{angle}$，反投影中的$N_{angle}$和$N_z$)，并使用11 GiB设备RAM作为GPU的可用内存，使用N³大小的图像和N²像素N个投影，那么GPU可以在投影耗尽内存之前处理N = 17000，而反投影则可以处理N = 8500。通过放宽限制并在每个GPU中要求单个图像切片和投影，限制可以增加到N = 27000。第一个限制将需要超过2TiB的CPU RAM来存储图像，而更宽松的GPU限制将需要74 TiB的CPU RAM。我们的结论是，任何这些大小的图像都不太可能在不久的将来使用，并且该代码可以安全地声称它可以在GPU上重建任意大的图像，即使存在上限。

推论是，在任何实际应用中，可能的图像大小的上限将由CPU RAM设置。除了将投影和图像完全存储在CPU RAM上的限制外，大多数迭代算法还需要相同大小的辅助变量。因此，所选的算法将设置现实上界。使用ROM作为交换内存或甚至只是写入和读取文件可能会增加映像大小超过可用的RAM，但会导致严重的性能瓶颈，因为辅助变量的持续读写ROM。我们决定不明确实现这一点。

在任何情况下，如果对预期运行的图像/投影大小添加约束，则有可能加速TIGRE中可用的一般断层扫描代码。如图9所示，不同的问题规模需要不同的优化策略。然而，我们提出了一个代码，在任何大小和数量的GPU的通用设置工作良好。对于具有特定架构和特定问题大小的给定GPU集，加速代码留给用户。

在多GPU方面，所提出的分割方法显示了单节点机器如何与多个GPU一起工作，然而，它仅限于这种计算机拓扑结构。幸运的是，该策略本身可以与多个节点的更高级别并行化相结合，类似于ASTRA工具箱所提供的[18]。使用MPI协议，代码可以适应多个节点拓扑，例如在HPC中。

图9中的结果显示了进一步可能的加速。对于512的大小，内存操作在投影中占用了大量的总时间，并且是反投影中的主要操作。然而，由于TIGRE的模块化设计，大多数这些操作都只存在。如果算法是直接在C++ /CUDA上编写的，那么这些操作在整个算法中只需要一次，而不是在每个投影和反投影调用中使用。因此，这样的实现将进一步加快算法的速度。这意味着所提供的代码可以用于时间紧迫的医疗案例，因为它可以实现重建时间短于一分钟。

计算结果表明，无论使用何种GPU，都可以成功地重建非常大的图像，并且可以非常快速地重建当前大小的图像。对于医学图像大小(通常小于512³)，在正确选择的硬件上，每次迭代不到1秒即可实现迭代重建。作为一个粗略的对比，TIGRE原文[26]表明，使用15次CGLS迭代可以在4min41s内完成512³张医学图像的重建。通过提出的实现，在单个GTX 1080 Ti上可以在1min01s内解决相同的问题。随着世界上图像尺寸的增加和CT机的使用范围的扩大，快速而不受资源限制的CT重建实现可以在医疗领域产生重大影响，允许任何拥有GPU的人使用迭代算法进行重建。

5. Conclusions

在这项工作中，我们提出了一种投影和反投影操作分割策略，该策略允许在单个节点的多个GPU上无缝分割计算，几乎没有开销。该方法对任意大的图像进行操作，消除了常见的GPU RAM大小限制，推高了单机CT可重构的最大尺寸的上界，并允许多台GPU机器显著加速。同样，它允许在具有小型GPU的机器中重建图像。此外，该方法允许更快地重建适合GPU的图像。这在医学CT中有明确的应用，在某些情况下重建时间是至关重要的。使用目前的实现，医疗尺寸图像的迭代重建可以在不到一分钟的时间内完成。