第一章 绪论
1.1 研究背景和意义
计算机断层成像技术(Computed Tomography, CT)产生于 20 世纪 70 年代，是一种由射线投影重建物体内部结构的成像技术[1]，具有扫描速度快、空间分辨率高等优点，是最重要的医学成像技术之一。
CT 图像重建算法分为解析法和迭代法。以滤波反投影(Filtered Back Projection, FBP)算法为代表的解析法是当前商用 CT 的主流重建方法。解析法重建速度快，但要求投影数据充分且完备。迭代法重建速度慢，但易于融入先验信息，形成基于压缩感知和稀疏先验的最优化重建方法。
CT 检查中过量的射线辐射会增加人体患癌风险，于是低剂量 CT 技术成为了当今 CT 算法研究的热点[2]。剂量的降低有两种方法，一种是只采集稀疏角度下的投影，一种是降低每个视角下投影所对应的射线剂量。前者所对应的重建可以被称为稀疏重建。然而解析法无法胜任稀疏重建，其重建的 CT 图像往往具有严重的条状伪影，无法被用于诊断评估。所以人们开始研究有望解决此问题的最优化重建方法。
2006 年，压缩感知[3–6]的提出为高精度迭代式稀疏重建提供了理论依据。从那时起，基于压缩感知的重建算法成为了稀疏重建的主流算法，大量优秀的算法被提出，极大地推动了稀疏重建的进展。全变分(Total Variation, TV)[7]最小算法是一种经典的压缩感知类重建算法，其简单有效，可以实现高精度稀疏重建。
TV 重建模型其实是一种最优化模型，包含一个数据保真项和一个 TV 正则项。可以认为 TV 算法是从满足数据约束的无穷多个解中选择出 TV 最小的解。数据保真项的作用是解成像模型所对应的线性方程组。该线性方程组的系数矩阵是投影矩阵，也称系统矩阵，在 CT 重建中往往异常巨大。
.........................
1.2 国内外研究现状
1970 年，代数重建算法的概念首次被 Gordon 等人引入图像重建领域中[8,9]。由于算法中投影和反投影巨大的计算量，导致代数重建算法的重建时间很长，如何能提高代数重建算法的重建效率，是当今 CT 图像重建研究的热点[10]。有很多提高代数重建效率的方法被提出。比如 Hudson 等人将投影数据分组为有序块，提出了有序子集排序方法[11]。Gordon 提出了一种无需正交便可重建的算法，投影数据以多级别访问的方式，将从上一级别中得到的投影数据充分利用在下一级别上[12]。蒋广胜等人通过将投影矩阵的方程均匀分配在曙光一号的四个处理机上，实现了代数重建法的多线程化，达到了很好的加速效率[13]。Mueller 提出了加权距离这一种新的投影排序方式，他优化了投影数据中的角距离，可对投影访问空间进行均匀采样，从而最小化投影序列中的相关性[14]。王小璞等人提出了一种基于分块迭代的快速代数重建算法，将图像逐级分块后，使用迭代算法将图像逼近，此方法需要实时计算投影系数矩阵的非零元[15]。除了改变投影数据的访问方式来提高算法的速度外，减少算法迭代的次数也可以达到此效果[16]。但是可以看到，无论哪种方式，都不可避免的需要投影矩阵这一数据，计算投影和反投影的工作量的本质并没有变化。
代数重建算法中，往往需要构建目标函数模型，通过特定的求解方程来得到迭代重建的公式[17]。目标函数模型一般包括保真项和正则项，因此人们利用此也可达到提高算法重建效率的目的。林少春等人将有序子集作为新的惩罚项引入加权最小二乘迭代算法，提高了算法的收敛速度[18]。Wang 等人提出了一种基于惩罚加权最小二乘(Penalized Weighted Least-Squares, PWLS)原理的迭代重建算法，设计了二次形式的各向异性惩罚项，有效保留了重建图像中的边缘信息[19]。1996 年，Rudin 等人根据偏微分方程的思想提出了全变分模型[20]，最早用于图像去噪领域。2006 年，Sidky 和Pan 等人在压缩感知理论基础上，提出了利用全变分最小化约束迭代重建的交替迭代最小化全变分(Projection Onto Convex Sets-Total Variation Minimization, POCS-TVM)算法[21]。2008 年，他们继续将扇形束 CT 重建扩展到三维锥形束 CT(Cone Beam Computed Tomography, CBCT)图像重建中，在 TV 最小化求解中引入可变步长，提出具有自适应的最速下降投影到凸集(Adaptive Steepest Descent Projection onto Convex Sets, ASD-POCS)算法，取得了不错的效果，可以有效抑制条形伪影，此方法避免了在傅里叶空间的运算[22]。
...........................
第二章 CT 成像原理与重建算法
2.1 CT 发展历史
计算机断层成像技术(Computed Tomography, CT)产生于 20 世纪 70 年代，是一种对通过射线对物体进行投影测量而获取物体内部组织结构的成像技术[1]，具有高时间分辨率、空间分辨率以及对比度分辨率的特点。其理论思想要追溯到 1917 年奥地利数学家 Radon 的贡献，他单单根据线形积分便确定了被积函数，其中积分指的是待重建图像的投影数据，被积函数为将要待重建的图像，实现了由投影重建图像这一数学设想，但是受限于当时的信息时代背景，这一理论一直被束之高阁，没有得到具体的应用。
1895 年，德国物理学家伦琴(Wilhelm Roentgen)因为发现了 X 射线而获得首届诺贝尔物理学奖。X 射线很快被应用到医学中，被用来探测人体疾病，成为医学诊断的常规手段。
然而经过 X 射线获得的人体影像是平面的，也就是说把立体的人体重叠在了一个平面上，这样获得到的信息是不完整的，无法获知人体内部的具体细节。后来人们发明了在不同角度拍摄的 X 射线断层摄影，在一定程度上弥补了原来技术的不足，但这仍然不能满足医学的需要。幸运的是，随着计算机技术的进步，出现了新的技术来弥补这一缺点[35]。
美国物理学家阿兰·麦克莱德·科马克(Allan Macleod Cormack)在研究放射性同位素 He6相关问题时意识到，人体内各部位组织是不均匀的，对于 X 射线的吸收或透过也是不均匀的，他把这个问题转化成了数学问题，并建立了数学模型。在 1957年和 1963 年，他制作了不同的样品作为模型来验证他的数学理论，并发表了相关的论文。不过他的论文起初并没有得到业界的注意，这主要是因为他的理论距离医学上的实践还有距离，同时当时的计算机技术也没有达到要求。
........................
2.2 CT 设备的基本组成
CT 设备主要由扫描系统、计算机系统和图像处理显示系统三部分组成。如图所示。扫描系统的作用是采集数据，X 射线管发出 X 光，穿过人体后被探测器吸收，光信号被转化为电信号进入模/数转化器(Analog/Digital Converter)，出来的数字信号被计算机系统存储。计算机系统的存储原则是将人体内存留的衰减系数重新排列成数字矩阵，数字矩阵被图像处理显示系统转化成不同灰度的体素，通过终端显示器显示最后的图像。
CT 系统不断更新换代，其获得图像的时间越来越短，且精度越来越高。按扫描方式划分，CT 技术先后经历了五代结构性能的发展和改变。
第一代 CT 机为平移-旋转扫描方式。最早一台于 1971 年安装在 winbledon 医院，这类扫描机多用于头部扫描，又被称为笔形扫描束装置。组成材料为一个 X 射线管和两到三个晶体探测器，X 射线管同探测器连接，X 射线管产生射线束穿过人体，同时探测器环绕人体做直线平移扫描运动。第一代扫描机虽然能够满足人体头部扫描，但是 X 射线利用率低且扫描时间长达 5 分钟，很难抑制图像伪影，故第二代扫描机出现。
第二代 CT 机的扫描方式依然是平移-旋转，但是它将 X 射线束换为了扇形线束，探测器数目也增加了数倍，扫描时间被缩短到 90 秒，但是图像的运动伪影仍不可避免，并且扇形线束在扩大扫描范围的同时也产生了更多的散射线，造成极大的浪费。
第三代 CT 机由美国通用电气公司(General Electric Company, GE)公司于 1975 年推出，它取消了平移运动而变为旋转-旋转扫描，扫描时间被缩短到 20 秒。同时三代采用较宽的扇形角，避免了 X 射线管产生的直线形焦点，这样也大大减少了伪影的产生[37]。目前的三代扫描已经可以将时间控制在 0.5 秒以内，在当今临床引用最为广泛。
................................
第三章 基于图像旋转的 TV 算法 .......................... 15
3.1 基于旋转的 ART 算法设计 .............................. 15
3.2 图像旋转中的插值算法 .......................... 16
第四章 基于图像旋转的非局部 TV 算法 ............................... 27
4.1 GPU 并行计算基础 .................................. 27
4.1.1 GPU 并行计算 .......................... 27
4.1.2 CUDA 架构 ................................ 28
第五章 总结与展望 ................................. 39
5.1 总结 ................................. 39
5.2 展望 ............................... 39
第四章 基于图像旋转的非局部 TV 算法
4.1 GPU 并行计算基础
在许多重要的科学和工程领域内，传统的计算方法往往无法满足实际的应用需求，因此出现了并行计算(Parallel Computing)。可以利用多种计算资源来处理计算问题，分为时间和空间两种，前者为流水线技术，后者则是运用多处理器来实现多指令流和数据流，其计算效率均高于串行计算(Serial Computing)，并行计算的多体系结构导致其实现依托于并行机和计算模型[40]。
4.1.1 GPU 并行计算
在图像处理领域，GPU 具有强大的并行计算能力，其浮点运算能力远远强于中央处理器(Central Processing Unit, CPU)[41,42]，也更加适应复杂的算法，因此人们常常通过 GPU 并行运算来实现算法的快速运行。不同于 CPU 中数量众多的控制单元和存储单元，GPU 的大部分空间被计算单元所占据，因此擅长大量并行线程和吞吐量密集的计算。
GPU 最初由 NVIDIA 公司提出，原型是三维图形流水线硬件，全称为 Graphics Processing Unit，由于可处理几乎所有的计算机图形的数据运算，故被称图形处理器。GPU 起初被局限于处理图形渲染等应用领域中，其图形渲染流水线的主要工作分为把三维坐标转化为二维坐标和把二维坐标转化为实际有颜色的像素这两个部分。在渲染更加复杂的三维场景图的过程中，研究人员偶然发现可以利用 GPU 的计算资源来进行一些并行计算的处理[43]，随着高性能计算设备的发展，GPU 已经由单纯的计算机图形处理转变成如今的并行计算工具。
..........................
5.1 总结
计算机断层成像技术作为现代医学中最为重要的诊疗手段，是医生精确诊断的得力工具，但其扫描过程中为患者带来健康风险的辐射剂量也不容被忽视，因此本文采用可以稀疏重建的迭代重建算法，提出了一种基于图像旋转的全变分迭代重建算法，主要研究内容如下：
1. 深入研究并实现了迭代重建算法，以及系统矩阵模型中的 Siddon 算法，为后续的算法研究做铺垫。
2. 提出了一种基于图像旋转的 TV 重建算法，该算法采用代数重建算法中的核心思想，在图像旋转的基础上，通过运用合理的插值算法，探索在不使用系统矩阵的情况下实现精确高效的迭代过程，并在迭代法的基础上，结合先验约束重建精确图像；通过仿真实验对提出的算法进行验证和评估，与其他经典 TV 算法比较了稀疏重建能力。实验表明，所设计的基于旋转的 TV 类算法在投影角度稀疏的情形下重建效果优于滤波反投影算法，与 ASD-POCS 算法相当。
3. 为了改善基于旋转的 TV 重建算法的精度，将算法的正则化部分改为 NLTV模型，同时使用 GPU 对 NLTV 模型权重矩阵的计算进行加速。设计了仿真实验对算法重建精度及加速效果进行了验证，实验表明，基于旋转的 NLTV 模型重建算法能有效改善 TV 模型的精度；通过 GPU 加速后的算法重建速度也有较大幅度提升。
参考文献（略）