并行医学图像重建：从图形处理器到网格

# 并行医学图像重建：从图形处理器到网格 在医学图像重建领域，并行计算技术的应用能够显著提升重建效率和质量。本文将详细介绍不同并行架构在医学图像重建中的应用，包括图形处理器（GPU）、Cell处理器等，并对它们的性能、编程舒适度、可访问性和成本效益进行比较。 ## 1. 图形处理器的并行化 ### 1.1 GPU架构与CUDA 现代图形处理器（GPU）可作为数学协处理器，为CPU增添计算能力。GPU是一种并行机器，由单指令多数据（SIMD）多处理器组成，其数量从1到32不等。SIMD多处理器的流处理器被称为着色器单元。GPU拥有自己的快速内存，容量可达4GB。在主板上，可同时安装1到4个GPU作为协处理器。例如，实验中使用的NVIDIA GeForce 8800 GTX，提供768MB设备内存，有16个多处理器，每个多处理器包含8个着色器单元。 NVIDIA通过CUDA（Compute Unified Device Architecture）为GPU提供了一个编程接口，将线程编程概念引入到C语言中。一组执行相同代码片段（即内核程序）的线程在一个多处理器上运行。该组中的每个线程在GPU的一个着色器单元上运行，每个单元在不同的数据元素上执行内核。应用程序的所有线程块由调度器分配到各个多处理器上。GPU的设备内存由所有线程共享。 ### 1.2 List - Mode OSEM CUDA实现 在GPU上计算一个子集的步骤如下： 1. CPU读取子集的事件，并将其复制到GPU设备内存中。 2. 每个线程计算$c_l$的部分和，并直接将其添加到设备内存中。每个线程处理的事件数量需综合考虑：一方面，应尽可能启动更多线程以有效隐藏内存延迟；另一方面，每个线程需在设备内存中保存部分结果，若为每个事件启动一个线程，会占用过多内存。因此，启动的线程数量应使所有部分结果仍能装入设备内存。 3. 每个线程计算$f_{l + 1}=f_lc_l$的一个体素值。 4. 将$f_{l + 1}$复制回CPU。 在计算$c_l$时，线程会像在共享内存实现中一样，直接写入共享向量$c_l$。为避免竞争条件，本应使用互斥锁保护$c_l$，但CUDA缺乏必要机制，只有原子整数操作。基于以下考虑，决定在GPU实现中允许竞争条件： - 当两个线程同时向一个体素添加一个浮点数时，最坏情况下一个线程会覆盖另一个线程的结果，即结果会稍有低估。 - 图像尺寸（如实验中的$150×150×280 = 6300000$）相较于并行写入线程的数量（如实验中的128）较大，因此竞争条件和错误体素的数量相对较少。实验估计，写入$c_l$时只有约0.04%会发生竞争条件。 - 所有体素的最大相对误差（由竞争条件和单精度浮点值的精度损失引起）小于1%，对重建图像无视觉影响。 - 核医学诊所中大多数小鼠和大鼠扫描的目的是确定特定器官（如肝脏）是否摄取了放射性物质。对于这些实验，重要的是重建图像的高精度，而不是精确的定量结果。对于定量实验，可使用线程安全的共享内存或Cell处理器进行重建。 当同时使用两个GPU时，有两个独立的设备内存。计算过程与上述相同，每个GPU在正向投影（步骤2）时计算一半的事件，在计算$f_{l + 1}$（步骤3）时计算一半的子图像。所有正向投影完成后，需对两个设备内存中的$c_l$进行求和。 ## 2. Cell处理器的并行化 ### 2.1 Cell架构及其编程 Cell宽带引擎是由索尼电脑娱乐公司、东芝公司和IBM公司联合开发的多处理器。它由一个PowerPC处理器元素（PPE）和八个称为协同处理器元素（SPE）的处理核心组成。通信通过元素互连总线（EIB）进行，包括PPE与SPE之间的通信、对共享内存（主存储）的访问以及I/O操作。PPE作为SPE的控制器，负责分配计算工作负载和处理操作系统任务。PPE由PowerPC处理器单元（PPU）和一个缓存组成。SPE通常负责处理程序的计算工作负载。每个SPE由一个运行频率为3.2 GHz的协同处理器单元（SPU）和一个256KB的本地存储（LS）组成，可通过DMA传输从主存储传输数据。DMA传输大小范围从128字节到16KB。SPU的本地存储数据访问速率为51GB/秒。 IBM为Cell处理器提供了软件开发工具包（SDK），其中包含用于PPU和SPU的GCC C/C++语言编译器。此外，SDK还包括一组Cell处理器C/C++库，为PPE端访问SPE提供应用程序编程接口（API）。在SPE端，库提供了用于在主存储和本地存储之间传输数据的DMA命令，包括提供互斥以避免竞争条件的原子命令。 ### 2.2 List - Mode OSEM Cell实现 在p个SPE上计算示例算法的一个子迭代步骤如下： 1. PPE读取子集的事件，并将其存储在主存储中。然后，PPE向每个SPE发送消息以启动计算。 2. 每个线程计算$c_l$的部分和，并直接将其添加到设备内存中。由于所有线程同时写入共享的$c_l$，使用原子操作。 3. 将重建图像划分为子图像$f^j$。每个SPE在其子图像上计算$f_{l + 1}^j=f_l^jc_l^j$。 在正向投影（步骤2）时，程序员需组织将$f_l$和$c_l$所需的体素从主存储传输到SPE的本地存储。由于DMA最小传输大小为128字节，计算$c_{l,j}$和$c_l + c_{l,j}$时，每个$f_l$的体素需传输128字节而非4字节。由于路径几乎都会穿过3D图像的多个y和z平面，传输的$f_l$的大部分额外体素在后续计算中无法使用。使用128字节的最小传输大小，每次DMA传输平均只有32个传输体素中的1.6个（即5%）被使用。 当使用两个Cell处理器（即两个PPE和总共16个SPE）时，有两个主存储，每个PPE仅与其SPE通信。两个主存储之间的通信和SPE管理对程序员是透明的，程序员可像只有一个带有16个SPE的PPE