步骤（2)，利用Grangeat公式将已知锥形束投影数据转化为已知三维Radon数据 关于p的一阶导数。Grangeat公式为：
[0084] 其中，幻足已知三维Radon数据关于P的一阶导数；如附图2、3、4所示，C〇si3 = S0/SCD，S0是放射源S与坐标原点0之间的距离，CD是放射线与探测器平面的交 点，则SCD是放射源S与交点C D之间的距离；P为坐标原点0与特征点C之间的距离；舜为 从坐标原点0指向特征点C的单位向量；;^/(S(p幻,?)是对步骤（1)中获取的不完备锥形 束投影数据进行加权计算后的投影数据值；所用的权值可优选为SO/SA，SA是 放射源S与直线t上点的距离；直线t与线段0CD垂直，且坐标原点0到直线t的垂直距离 为s，p，q分别为探测器平面的横坐标轴和纵坐标轴，α为线段0CD与探测器平面横坐标轴 P的夹角；
[0085] 然后对已知三维Radon数据关于P的一阶导数/?'/'(ρ??)进行积分即可获得已知三 维Radon数据&瓦)〇:
[0086] 步骤（3)，计算已知三维Radon数据财(两)的投影几何矩变换公式为：
[0088] 其中，为在方向S上的P阶投影几何矩变换。
[0089] 通过式（2)计算阶数P分别为[0，1，2，···，Μ]的一系列投影几何矩变换，Μ为所使 用的几何矩的最大阶数。在本实例中Μ =20。然后，将这些投影几何矩变换值按一定顺序 组成向量Φμ<Χ):
[0091] 步骤（4)，针对已扫描角度范围，计算系数矩阵
[0094] 步骤（5)，已知三维Radon数据的投影几何矩变换与几何图像矩之间关系的矩阵 形式为：
[0096] 则利用公式（5)可得几何图像矩向量为：
[0098] 几何图像矩向量ΨΜ由一系列阶数不超过Μ的几何图像矩G_组成，如下所示：
[0101] 步骤（6)，针对未扫描角度范围，计算系数矩阵￡1?〇〇:
[0103] 其中，_为未扫描角度范围内的向量，BP
[0104] 步骤（7)，未知Radon数据的投影几何矩变换与几何图像矩之间关系的矩阵形式 为：
[0106] 因为步骤（5)计算出的几何图像矩ΨΜ与拜、皮均无关，所以可将步骤（5)获得的 几何图像矩ΨΜ以及步骤（6)获得的系数矩阵:!QM(rY)代入公式（9)，即可计算出未知Radon 数据的投影几何矩变换$_〇〇, 是一系列阶数不大于Μ的投影几何矩组成 的向量（〇彡P彡M)。
[0107] 步骤（8)，投影几何矩逆变换公式如下：
[0109] 则利用公式（10)和步骤（7)获得的Zp〇〇即可求得未知Radon数据
[0110] 步骤（9)，将估计出的未知Radon数据按照扫描角度顺序和已知Radon数据合并， 即可获得补全的Radon数据。
[0111] 步骤（10)，三维Radon逆变换公式如下：
[0113] 利用三维Radon逆变换即可从补全的Radon数据中重建出CT图像/(幻
[0114] 为了对比本发明所示方法的效果，使用传统的滤波反投影重建算法对不完备锥形 束投影数据进行重建，得到的重建图像如图5所示；而使用本发明方法获得的重建图像如 图6所示。比较图5和图6的重建结果，可以看出本发明方法所得结果中伪影相对较少，而 且头模型的轮廓更加完整和清晰。因此可以认为本发明方法能够在部分投影数据缺失的条 件重建出伪影更少、图像质量更好的结果。
[0115] 以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是 按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明的保护范 围之内。
【主权项】
1. 基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在于，包括以下步骤： (1) 获取有限角度扫描条件下的不完备的锥形束投影数据； (2) 利用Grangeat公式将步骤（1)获得的锥形束投影数据转化为已知三维Radon数 据； (3) 对步骤（2)获得的已知三维Radon数据进行投影几何矩变换； (4) 建立步骤（3)获得的已知三维Radon数据的投影几何矩变换与几何图像矩之间的 关系式； (5) 利用步骤（4)建立的关系式，从步骤（3)获得的已知三维Radon数据的投影几何矩 变换中计算出几何图像矩； (6) 建立未知三维Radon数据与几何图像矩之间的关系式； (7) 根据步骤（6)建立的关系式，从步骤（5)获得的几何图像矩中估计出未知三维 Radon数据的投影几何矩变换； (8) 利用投影几何矩逆变换从步骤（7)获得的未知三维Radon数据的投影几何矩变换 中求出未知三维Radon数据； (9) 将步骤⑵得到的已知三维Radon数据和步骤⑶得到的未知三维Radon数据进 行数据拼合，获得补全的三维Radon数据； (10) 通过三维Radon逆变换，根据步骤（9)获得的补全的三维Radon数据重建出CT图 像。2. 如权利要求1所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于：在步骤（1)中，所述有限角度扫描条件下是指在[<p,范围内旋转扫描，其中， 0 ()<φ<360°，3. 如权利要求2所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于，步骤（2)的具体过程如下： 首先利用Grangeat公式将锥形束投影数据转化为已知三维Radon数据关于P的一阶 导数：式（1)中，是已知三维Radon数据关于P的一阶导数，cos β = S0/SCD，Cd是 放射线与探测器平面的交点，则SCd是放射源S与交点C D之间的距离，SO是放射源S与坐 标原点0之间的距离，P为坐标原点0与特征点C之间的距离，S为从坐标原点0指向特征 点C的单位向量，是对步骤（1)中获取的不完备锥形束投影数据进行加权计 算后的投影数据值，直线t与线段OCd垂直，且坐标原点0到直线t的垂直距离为s，p，q分 别为探测器平面的横坐标轴和纵坐标轴，α为线段OCd与探测器平面横坐标轴p的夹角； 然后对已知三维Radon数据关于P的一阶导数进行积分，即获得已知三维 Radon数据，所述三维Radon数据定义为：式⑵中，/(幻是三维图像f在添处的灰度值，Θ 是单位向量海__与坐标轴Z的夹角。4. 如权利要求3所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于：将SO/SA作为权值，对步骤（1)中获取的不完备锥形束投影数据进行加权计算后得到投 影数据值其中，SA是放射源S与直线t上任意点A的距离。5. 如权利要求3所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于，所述步骤（3)中的投影几何矩变换的定义为：式（3)中，Lp 〇i)为在方向i上的p阶投影几何矩变换数据。6. 如权利要求5所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于，所述步骤（4)中的几何图像矩的定义为：式（4)中，f (X，y, z)是三维图像f在点（X，y, z)处的灰度值。 则已知三维Radon数据的p阶投影几何矩变换与几何图像矩之间的关系式为：7. 如权利要求6所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于，所述步骤（5)的具体过程如下： 首先将式（5)改写成如式（7)所示的矩阵形式：M为所用几何矩的最大阶数； 然后对式（7)采用矩阵除法，获得几何图像矩向量：8. 如权利要求6或7所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征 在于，所述步骤(6)中的建立几何图像矩与未知三维Radon数据的投影几何矩变换之间的 关系式为：式（12)中，IpOO是未知Radon数据的投影几何矩变换，为未扫描角度范围内的向 量，即9. 如权利要求8所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于，所述步骤（7)的具体过程为： 首先将式（12)改写成如式（13)所示的矩阵形式：然后将步骤（5)获得的几何图像矩向量ΦΜ代入式（13)，计算出未知三维Radon数据 的投影几何矩变換10. 如权利要求9所述基于几何图像矩的有限角度锥形束CT图像重建方法，其特征在 于，所述步骤（8)中投影几何矩逆变换的公式为：其中，是未知三维Radon数据； 所述步骤（10)中三维Radon逆变换的公式为：I. 其中，八幻是重建的CT图像。
【专利摘要】本发明公开了基于几何图像矩的有限角度锥形束CT图像重建方法，步骤依次为，将获得的锥形束投影数据转化为已知Radon数据、获得已知Radon数据的投影几何矩变换、根据已知Radon数据的投影几何矩变换中计算出几何图像矩、根据几何图像矩估计出未知Radon数据的投影几何矩变换、经逆变换求出未知Radon数据、将步已知Radon数据和未知Radon数据进行数据拼合，获得补全的三维Radon数据并重建CT图像。本发明能够在减小扫描范围的条件下重建出符合临床诊断要求、高质量的锥形束CT图像。
【IPC分类】G06T15/06, G06T11/00
【公开号】CN105354868
【申请号】CN201510644674
【发明人】戴修斌, 刘天亮, 胡栋, 晏善成, 石丹丹, 邓黄建
【申请人】南京邮电大学
【公开日】2016年2月24日
【申请日】2015年10月8日