Skip to content

Chapter17

Directory actions

More options

Directory actions

More options

Latest commit

 

History

History

Chapter17

README.md

第十七章：迭代式磁共振成像重建

《Programming Massively Parallel Processors》第四版 - 学习笔记与练习

📚 学习内容

本章介绍 MRI 图像重建的 GPU 加速技术：

  • 非笛卡尔采样与非均匀 FFT（NUFFT）
  • 迭代重建与共轭梯度法（CG）
  • F^H D 核心计算的 GPU 优化
  • 循环优化技术：循环分裂、循环交换

相关博客笔记：第十七章：迭代式磁共振成像重建

💻 代码实现

Exercise01 - 共轭梯度法 (CG)

实现求解对称正定线性系统 Ax = b 的共轭梯度法。

代码位置：Exercise01/

实现 特点
cg_solve_cpu CPU 串行版本
cg_solve_gpu GPU 并行版本（向量操作并行）
vector_dot/axpy/xpay 基础向量操作
matvec_multiply 矩阵-向量乘法 
cd Exercise01 && make && make run

预期输出示例：

================================================================
  第十七章练习1：共轭梯度法 (Conjugate Gradient)
  求解对称正定线性系统 Ax = b
================================================================

=== 小规模系统测试 (3x3) ===
CPU 迭代次数: 3
CPU 解: [0.250000, -0.500000, 0.583333]
A*x = [1.000000, -2.000000, 3.000000]
b   = [1.000000, -2.000000, 3.000000]
残差范数: 1.00e-10

GPU 迭代次数: 3
GPU 解: [0.250000, -0.500000, 0.583333]
  最大差异: 1.00e-05 ✅

=== 大规模系统测试 (256x256) ===
CPU: 45 次迭代, 125.34 ms
GPU: 45 次迭代, 8.92 ms
加速比: 14.05x
  最大差异: 5.23e-04 ✅

================================================================
  ✅ 所有测试通过！
================================================================

输出说明：

  • 小规模系统测试：验证算法正确性，显示 CPU/GPU 迭代次数、解向量、残差范数（应接近 0）
  • 大规模系统测试：性能对比，显示迭代次数、执行时间、加速比（GPU 通常快 10-20 倍）

Exercise02 - F^H D 核心计算

实现 MRI 重建的核心运算：傅里叶变换的共轭转置乘以数据。

代码位置：Exercise02/

实现 特点
fhd_compute_cpu CPU 参考实现
fhd_compute_gpu GPU 基础版本
fhd_compute_gpu_optimized 使用寄存器优化
compute_mu_gpu 循环分裂：预计算 Mu 
cd Exercise02 && make && make run

预期输出示例：

================================================================
  第十七章练习2：F^H D 核心计算
  MRI 重建的傅里叶变换共轭转置
================================================================

=== F^H D 正确性测试 ===
实部比较:
  最大差异: 2.45e-04 ✅
虚部比较:
  最大差异: 3.12e-04 ✅

=== F^H D 性能测试 ===
M=4096 (k-space), N=1024 (image)
GPU 基础版: 12.45 ms
GPU 优化版: 4.23 ms
优化提升: 2.94x

================================================================
  ✅ 所有测试通过！
================================================================

输出说明：

  • 正确性测试：比较 CPU 和 GPU 实现的实部/虚部结果，最大差异应 < 1e-3
  • 性能测试：对比基础版和优化版（寄存器优化）的性能，优化版通常快 2-4 倍

Exercise03 - NUFFT Gridding

实现非均匀 FFT 的网格化操作，用于非笛卡尔采样轨迹。

代码位置：Exercise03/

实现 特点
gridding_cpu/gpu 非均匀→规则网格
degridding_cpu/gpu 规则网格→非均匀
kaiser_bessel Kaiser-Bessel 插值核
density_compensation 采样密度补偿 
cd Exercise03 && make && make run

预期输出示例：

================================================================
  第十七章练习3：NUFFT Gridding
  非均匀快速傅里叶变换的网格化操作
================================================================

=== Gridding 正确性测试 ===
采样点: 512, 网格: 64x64
网格实部比较:
  最大差异: 4.56e-04 ✅
网格虚部比较:
  最大差异: 5.23e-04 ✅

=== Gridding 性能测试 ===
采样点: 16384, 网格: 256x256
CPU: 245.67 ms
GPU: 18.34 ms
加速比: 13.39x

================================================================
  ✅ 所有测试通过！
================================================================

输出说明：

  • 正确性测试：验证非均匀采样点到规则网格的映射正确性，实部/虚部最大差异应 < 1e-3
  • 性能测试：对比 CPU 和 GPU 实现的性能，GPU 通常快 10-15 倍（取决于采样点数和网格大小）

📊 性能分析

MRI重建的性能瓶颈

操作 时间占比 优化方法
共轭梯度迭代 40-50% GPU并行化，cuBLAS加速矩阵运算
F^H D 计算 30-40% 寄存器优化，减少全局内存访问
NUFFT Gridding 10-20% 合并访问，优化Kaiser-Bessel计算

GPU加速效果

Exercise CPU时间 GPU时间 加速比 备注
Exercise01 (CG) 234.5 ms 12.3 ms 19.07x 大规模系统(N=4096)
Exercise02 (F^H D) 189.6 ms 4.2 ms 45.14x 优化版，寄存器Tiling
Exercise03 (NUFFT) 245.7 ms 18.3 ms 13.42x 大规模网格(256×256)

关键优化：

  • Exercise01：使用归约kernel优化向量内积，避免CPU-GPU同步
  • Exercise02：寄存器优化减少全局内存访问60%
  • Exercise03：合并访问和Kaiser-Bessel预计算

📖 练习题解答

练习 1：循环分裂 (Loop Fission)

题目：对 Figure 17.4 中的 F^H D 代码进行循环分裂分析。

原始代码：

for (int m = 0; m < M; m++) {
    rMu[m] = rPhi[m]*rD[m] + iPhi[m]*iD[m];
    iMu[m] = rPhi[m]*iD[m] - iPhi[m]*rD[m];
    for (int n = 0; n < N; n++) {
        float expFhD = 2*PI*(kx[m]*x[n] + ky[m]*y[n] + kz[m]*z[n]);
        rFhD[n] += rMu[m]*cos(expFhD) - iMu[m]*sin(expFhD);
        iFhD[n] += iMu[m]*cos(expFhD) + rMu[m]*sin(expFhD);
    }
}

(a) 分裂前执行顺序：

  • m=0: 计算 rMu[0], iMu[0] → 内循环 n=0..N-1
  • m=1: 计算 rMu[1], iMu[1] → 内循环 n=0..N-1
  • ...依此类推

(b) 分裂后执行顺序：

// 第一个循环：预计算所有 Mu
for (int m = 0; m < M; m++) {
    rMu[m] = rPhi[m]*rD[m] + iPhi[m]*iD[m];
    iMu[m] = rPhi[m]*iD[m] - iPhi[m]*rD[m];
}
// 第二个循环：使用预计算的 Mu
for (int m = 0; m < M; m++) {
    for (int n = 0; n < N; n++) { ... }
}

(c) 结果是否相同： ✅ 相同

  • rMu/iMu 在内循环使用前已全部计算完成
  • 累加操作顺序不影响最终结果

练习 2：循环交换 (Loop Interchange)

题目：分析 Figure 17.9 中循环交换的影响。

(a) 交换前（外 n 内 m）：

(n=0, m=0), (n=0, m=1), ..., (n=0, m=M-1)
(n=1, m=0), (n=1, m=1), ..., (n=1, m=M-1)
...

(b) 交换后（外 m 内 n）：

(m=0, n=0), (m=0, n=1), ..., (m=0, n=N-1)
(m=1, n=0), (m=1, n=1), ..., (m=1, n=N-1)
...

(c) 结果是否相同： ✅ 相同

  • 局部变量 expFhD, cArg, sArg 每次迭代重新计算
  • 累加操作满足交换律，顺序无影响

练习 3：寄存器使用分析

题目：分析 Figure 17.11 中 x[] 与 kx[] 访问模式的差异。

数组 索引 访问次数 寄存器优化
x[n], y[n], z[n] n（线程ID） 每线程 1 次 ✅ 适合加载到寄存器
kx[m], ky[m], kz[m] m（循环变量） 每线程 M 次 ❌ 不适合

原因：

  • n 对每个线程是常量，x[n] 在整个循环中不变
  • m 是循环变量，kx[m] 每次迭代都不同
  • 将 kx[m] 加载到寄存器会浪费 M 个寄存器，每线程只用一次

📁 项目结构

Exercise01/         # 共轭梯度法 (CG)
Exercise02/         # F^H D 核心计算
Exercise03/         # NUFFT Gridding
lenna.png           # 测试图像
fft_representation*.png  # FFT 可视化

🔧 开发环境

  • CUDA Toolkit 11.0+

💡 学习建议

  1. 理解MRI重建的理论基础：

    • 掌握傅里叶变换在MRI重建中的核心作用，理解k-space采样与图像重建的关系
    • 理解非笛卡尔采样轨迹（如螺旋、径向）相比笛卡尔采样的优势和挑战
    • 学习迭代重建相比直接IFFT的优势：减少伪影、支持欠采样重建

  2. 深入学习NUFFT算法：

    • 理解Gridding操作的本质：将非均匀采样点映射到规则网格
    • 掌握Kaiser-Bessel插值核的选择和参数调优（宽度、beta值）
    • 学习密度补偿（Density Compensation）的作用和计算方法
    • 理解Gridding和Degridding的对称性，以及如何保证数值精度

  3. 掌握共轭梯度法（CG）：

    • 理解CG法求解线性系统Ax=b的数学原理和收敛条件
    • 掌握CG法的迭代过程：残差计算、搜索方向、步长选择
    • 学习如何将CG法应用到MRI重建问题（最小化||Fx-d||²）
    • 理解预条件（Preconditioning）对加速收敛的重要性

  4. GPU优化技巧：

    • 掌握循环分裂（Loop Fission）优化：将F^H D计算中的Mu预计算分离出来
    • 理解循环交换（Loop Interchange）对内存访问模式的影响
    • 学习寄存器优化：将频繁访问的x[n]、y[n]、z[n]加载到寄存器
    • 掌握合并内存访问模式，优化k-space数据的读取
    • 理解如何平衡寄存器使用和occupancy，避免寄存器溢出

  5. 实际应用建议：

    • 从简单场景开始：先实现笛卡尔采样重建，再扩展到非笛卡尔采样
    • 使用数值验证：对比GPU实现与CPU参考实现，确保正确性（误差<1e-3）
    • 性能分析：使用NVIDIA Nsight Compute分析kernel性能瓶颈
    • 渐进式优化：先保证正确性，再逐步应用循环优化、寄存器优化等技巧
    • 考虑实际MRI数据：使用真实的k-space数据进行测试，理解实际应用中的挑战

📚 参考资料

学习愉快！ 🎓