GPU Gems 2 · Chapter 48. Medical Image Reconstruction with the FFT

本章讨论的不是通用的渲染技巧清单，而是一个很具体的 GPU 数值计算问题：如何把医学成像中的 Fourier-domain reconstruction 变成可并行、可检查的纹理数据流。NVIDIA 原章由 Thilaka Sumanaweera 和 Donald Liu（Siemens Medical Solutions USA）撰写，分别以 MRI 与超声说明 GPU FFT 的价值。

1. 为什么医学图像需要 Fourier-domain reconstruction

MRI 与超声都先得到仪器采集的信号，再从信号恢复空间图像。原始信号往往更适合在频率域组织：MRI 采集的是与质子密度 Fourier transform 相关的 k-space 数据；超声 PPI 则从不同换能器位置记录回波。重建的计算量集中在二维 FFT、复数乘法、频率坐标变换与逆变换。

Fourier transform 让我们从另一套坐标观察同一份信息。它不等于“压缩成一张频谱图”：幅度和相位都重要，医学重建通常还要保留复数值，最后用 inverse transform 回到图像域。

这里有一个工程判断：如果只需要解释低频趋势，可以观察幅度；如果要重建边缘、位置和相位，必须保持复数分量以及采样坐标。GPU 加速的是大量规则的复数运算，不会自动修复缺失采样或不正确的物理模型。

2. 从 DFT 到 FFT：相同结果，不同计算组织

discrete Fourier transform（DFT）可以直接写成“每个输出频率检查所有输入样本”的双重循环。对 N 个样本，它需要约 N² 次复数级别的工作；在二维图像上，两个维度会让代价更明显。

for k in outputFrequencies:
  spectrum[k] = 0
  for n in inputSamples:
    angle = -2π * k * n / N
    spectrum[k] += input[n] * exp(i * angle)

FFT 的结果与 DFT 相同，改变的是计算顺序。它把偶数和奇数索引拆开，复用子问题，再把局部结果通过 twiddle weight 合并。对适合的 N，复杂度从 O(N²) 降为 O(N log₂N)。

butterfly 是这个分解的局部单元。对 radix-2 FFT，长度通常取 2 的幂；以 8 点输入为例，需要先做 bit-reversal，再执行 2 点、4 点、8 点三层 butterfly。复杂度下降来自重复子问题被复用，不是来自降低结果精度。

3. 把 FFT stage 映射到 GPU multipass

二维 FFT 可以先对每一列做一维 FFT，再对每一行做一维 FFT。GPU Gems 2 的实现思路是把每个一维 FFT stage 写成 fragment pass：fragment 根据纹理坐标找到两个输入，乘以 stage 对应的复数权重，再写入两个输出位置。

关键不是“把循环翻译成 shader”这么简单，而是保留 stage 的数据依赖。每个 pass 都应该明确：当前读哪一个 source、结果写哪一个 target、下一层何时交换、复数的实部和虚部放在哪些通道。原章还指出，可以把第一轮 twiddle weight 为 1 的计算和 input scrambler 合并，并用同一个 pbuffer 的多个 draw buffers 承接中间结果，避免 CPU 往返。

bindDrawBuffer(source);
bindDrawBuffer(target);
setUniform("stage", stage);
setUniform("size", transformSize);
drawFullscreenQuad();
swap(source, target);

这段伪代码表达的是资源契约，不是某个现代 API 的可复制实现。验证时可以保存每层的中间纹理，与 CPU 参考 FFT 比较；如果 stage 1 就不一致，优先查 bit-reversal 和 source/target，而不是先调线程数。

4. MRI：k-space、二维 IFFT 与空间图像

k-space 可以理解为 MRI 图像的空间频率坐标系。经过射频信号解调后，采集结果与目标区域质子密度的 Fourier transform 相关；扫描过程沿 kx、ky 逐线穿过频率空间，完整或部分采样后再用二维 inverse FFT 重建。

实践中需要明确四件事：采样顺序、复数共轭与相位约定、k-space 是否中心化，以及 zero-padding 或缺失采样如何处理。GPU FFT 只负责规则的变换内核；采样校正、窗函数、噪声处理和临床显示仍是上游或下游的责任。

5. 超声 PPI：二维 FFT 后的频率重映射

在 pulse plane-wave imaging（PPI）中，换能器阵列可以同时发射平面波，记录空间位置 x 与时间 t 上的回波 s(x,t)。对这个二维数据做 FFT 得到 S(fx,ft) 后，传播关系把时间频率重新解释为深度频率，再把结果送入二维 inverse FFT。

frequency remapping 是超声路径与 MRI 路径最容易混淆的区别。它通常需要插值：某个 (fx, ft) 不一定正好落在目标 (fx, fz) 网格上，因此要处理边界、空洞、幅度权重与相位。remap 前后的坐标轴必须写在调试图上，不能只凭最终图像猜测。

如果直接把 S(fx,ft) 当作 S(fx,fz) 做 IFFT，算法仍可能输出一张“像图像”的结果，但深度位置和结构会错误。正确的诊断顺序是：检查频率坐标，再检查插值覆盖率，接着检查复数相位与 IFFT 缩放，最后才比较 GPU 与 CPU 性能。

6. 先预测，再观察 FFT 的中间层

先预测：选择 MRI 信号，在 stage 0 时可以看到 8 个原始样本；推进到 bit-reversal 后，数值顺序会改变但总能量不应凭空消失；推进到最后一层后，低频与周期结构会集中到特定频率分量。切换超声信号，预测峰值会如何移动，再用实部/幅度切换检查“相位信息不能只看柱高”这一点。

这个实验实际执行 8 点 radix-2 复数 FFT 的重排、twiddle weight 乘法、加法和减法。它把 stage 暴露出来，是为了让 source/target 的边界和复杂度推理可观察；它不是医学诊断模拟，也不包含真实扫描噪声、采样轨迹或临床校准。

7. 正确性与性能边界

工程上可以用以下检查表收敛 GPU FFT：

1. 数值参考：用同一输入、同一正负号和同一归一化约定，与 CPU DFT 或可信 FFT 库逐 stage 比较。
2. 布局参考：确认行 FFT 和列 FFT 的纹理坐标映射；转置、stride、复数通道和 bit-reversal 任何一项错位都会污染全图。
3. 资源参考：每层使用独立 source/target，避免读写 hazard；二维处理完成后再进入 remapping 或逆变换。
4. 误差参考：单独记录最大绝对误差、相位误差、插值空洞和边界能量，不要用一张颜色图代替数值断言。
5. 性能参考：GPU 的收益来自并行 butterfly、纹理带宽和减少 CPU-GPU 往返；原章的硬件性能数字属于历史实验，不应直接当作今天的基准承诺。

for each row:    FFT(row, source, target)
transpose or map coordinates
for each column: FFT(column, source, target)
if ultrasound:   remap(fx, ft -> fx, fz)
image = inverseFFT(source, target)

小结

• Fourier transform 改变信号的坐标表示；DFT 直接求和，FFT 通过 butterfly 复用子问题。
• GPU multipass FFT 必须显式处理 bit-reversal、stage 顺序、复数布局和 source/target 交换。
• MRI 以 k-space 为中间表示，二维 inverse FFT 把空间频率采样带回图像坐标。
• 超声 PPI 需要先对 s(x,t) 做二维 FFT，再进行 frequency remapping，最后 inverse FFT。
• GPU 加速数值内核，但不会替代采样校正、物理模型、插值质量和 CPU 参考验证。