基于几何图像矩的有限角度锥形束CT图像重建方法

摘要

本发明公开了基于几何图像矩的有限角度锥形束CT图像重建方法，步骤依次为，将获得的锥形束投影数据转化为已知Radon数据、获得已知Radon数据的投影几何矩变换、根据已知Radon数据的投影几何矩变换中计算出几何图像矩、根据几何图像矩估计出未知Radon数据的投影几何矩变换、经逆变换求出未知Radon数据、将步已知Radon数据和未知Radon数据进行数据拼合，获得补全的三维Radon数据并重建CT图像。本发明能够在减小扫描范围的条件下重建出符合临床诊断要求、高质量的锥形束CT图像。

说明书

技术领域

本发明涉及一种医学图像重建方法，特别涉及了一种CT图像重建方法。

背景技术

作为目前一种常规有效的临床医学诊断工具，X射线计算机断层成像技术(X-ray ComputerizedTomography,CT)为临床医生的诊断提供了丰富的人体器官组织信息。但是由 相关研究表明：一次完整的CT扫描通常伴随着较高程度的电离辐射，而高剂量电离辐射 可诱发人体新陈代谢异常乃至癌症、白血病等疾病。因此，如何在降低X射线使用剂量的 同时，保证重建图像质量满足临床诊断要求成为医学图像处理领域研究的重点。

临床上减少病患辐射量的重要方法之一就是减小CT扫描范围，即将探测器的旋转角 度范围限制在某个小于标准的区间内，从而在总体上大幅减少了患者所受X射线辐射量。 虽然限制CT设备扫描范围能够降低患者所受X射线辐射量，但同时会造成所获CT投影 数据部分缺失，即获得的是不完备投影数据，使重建CT图像质量明显下降，以至于无法 满足临床诊断的需要。

目前针对有限角度CT图像重建方法主要是基于统计模型的迭代重建方法。虽然该方 法在扫描角度有限的情况下仍能获得较好的重建结果，但是该方法需要对目标函数进行上 百次的迭代求解，重建CT图像所需计算时间大幅增加，远远超过经典的解析重建方法， 不能满足临床CT实时显像的要求。

经典的解析重建方法虽然重建速度快，但是在扫描角度有限的情况下会丢失图像原有 细节信息，从而导致重建图像中出现大量伪影和噪声，严重影响了对特征点的分辨。

发明内容

为了解决上述背景技术提出的技术问题，本发明旨在提供基于几何图像矩的有限角度 锥形束CT图像重建方法，能够在减小扫描范围---即投影数据不完备条件下重建出符合临 床诊断要求、高质量的锥形束CT图像。

为了实现上述技术目的，本发明的技术方案为：

基于几何图像矩的有限角度锥形束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)将步骤(2)得到的已知三维Radon数据和步骤(8)得到的未知三维Radon数据进 行数据拼合，获得补全的三维Radon数据；

(10)通过三维Radon逆变换，根据步骤(9)获得的补全的三维Radon数据重建出CT图 像。

进一步地，在步骤(1)中，所述有限角度扫描条件下是指在范围内旋转扫描， 其中，

进一步地，步骤(2)的具体过程如下：

首先利用Grangeat公式将锥形束投影数据转化为已知三维Radon数据关于ρ的一阶导数：

$R^{\prime }f\left(\rho \bar{n}\right)=\frac{1}{{\cos }^2\beta }\left(cos\alpha {\int }_{-\infty }^{\infty }\frac{\partial }{\partial p}{X}_{w}f\right(s\left(\rho \bar{n}\right),t)dt+\sin \alpha {\int }_{-\infty }^{\infty }\frac{\partial }{\partial q}{X}_{w}f(s\left(\rho \bar{n}\right),t\left)dt\right)---\left(1\right)$

式(1)中，是已知三维Radon数据关于ρ的一阶导数，cosβ＝SO/SC_D，C_D是放射线 与探测器平面的交点，则SC_D是放射源S与交点C_D之间的距离，SO是放射源S与坐标原点O之 间的距离，ρ为坐标原点O与特征点C之间的距离；为从坐标原点O指向特征点C的单位向 量；是对步骤(1)中获取的不完备锥形束投影数据进行加权计算后的投影数 据值，直线t与线段OC_D垂直，且坐标原点O到直线t的垂直距离为s，p,q分别为探测器平面 的横坐标轴和纵坐标轴，α为线段OC_D与探测器平面横坐标轴p的夹角；

然后对已知三维Radon数据关于ρ的一阶导数进行积分，即获得已知三维Radon数据， 所述三维Radon数据定义为：

$Rf\left(\rho \bar{n}\right)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f\left(\bar{x}\right)\delta (\bar{x}·\bar{n}-\rho )d\bar{x}---\left(2\right)$

式(2)中，是三维图像f在处的灰度值，θ是单位 向量与坐标轴z的夹角。

进一步地，将SO/SA作为权值，对步骤(1)中获取的不完备锥形束投影数据进行加权计算 后得到投影数据值其中，SA是放射源S与直线t上任意点A的距离。

进一步地，所述步骤(3)中的投影几何矩变换的定义为：

$L_p\left(\bar{n}\right)={\int }_{-\infty }^{\infty }\rho ·{\mathrm{Rf}}_{\bar{n}}\left(\rho \right)d\rho ---\left(3\right)$

式(3)中，为在方向上的p阶投影几何矩变换数据。

进一步地，所述步骤(4)中的几何图像矩的定义为：

$G_{ijk}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{x}^i{y}^j{z}^{k}f\left(x,y,z\right)dxdydz---\left(4\right)$

式(4)中，f(x,y,z)是三维图像f在点(x,y,z)处的灰度值。

则已知三维Radon数据的p阶投影几何矩变换与几何图像矩之间的关系式为：

$L_p\left(\bar{n}\right)={\Sigma }_{m=0}^p{\Sigma }_{r=0}^{p-m}{ϵ}_{mr}(p,\bar{n})G_{m,r,p-m-r}---\left(5\right)$

其中，

进一步地，所述步骤(5)的具体过程如下：

首先将式(5)改写成如式(7)所示的矩阵形式：

${\Phi }_M\left(\bar{n}\right)={\Omega }_M\left(\bar{n}\right){\Psi }_M---\left(7\right)$

其中，

${\Phi }_M\left(\bar{n}\right)={[L_0\left(\bar{n}\right),L_1\left(\bar{n}\right),...L_M\left(\bar{n}\right)]}^T---\left(8\right)$

${\Psi }_M={[G^{{\left(0\right)}^T},G^{{\left(1\right)}^T},...G^{{\left(M\right)}^T}]}^T$

$G^{\left(k\right)}={[G_{k,0,M-k},G_{k-1,1,M-k},...,G_{1,k-1M-k},G_{0,k,M-k}]}^T---\left(9\right)$

M为所用几何矩的最大阶数；

然后对式(7)采用矩阵除法，获得几何图像矩向量：

${\Psi }_M={\Phi }_M\left(\bar{n}\right)/{\Omega }_M\left(\bar{n}\right)---\left(11\right)$

进一步地，所述步骤(6)中的建立几何图像矩与未知三维Radon数据的投影几何矩变换之 间的关系式为：

${\tilde{L}}_p\left({\bar{n}}^{\prime }\right)={\Sigma }_{m=0}^p{\Sigma }_{r=0}^{p-m}{ϵ}_{mr}(p,{\bar{n}}^{\prime })G_{m,r,p-m-r}---\left(12\right)$

式(12)中，是未知Radon数据的投影几何矩变换，${\bar{n}}^{\prime }$为未扫描角度范围内的向量，即

进一步地，所述步骤(7)的具体过程为：

首先将式(12)改写成如式(13)所示的矩阵形式：

${\tilde{\Phi }}_M\left({\bar{n}}^{\prime }\right)={\Omega }_M\left({\bar{n}}^{\prime }\right){\Psi }_M---\left(13\right)$

其中，

${\tilde{\Phi }}_M\left({\bar{n}}^{\prime }\right)={[L_0\left({\bar{n}}^{\prime }\right),L_1\left({\bar{n}}^{\prime }\right),...L_M\left({\bar{n}}^{\prime }\right)]}^T---\left(14\right)$

${\Psi }_M={[G^{{\left(0\right)}^T},G^{{\left(1\right)}^T},...G^{{\left(M\right)}^T}]}^T$

$G^{\left(k\right)}={[G_{k,0,M-k},G_{k-1,1,M-k},...,G_{1,k-1,M-k},G_{0,k,M-k}]}^T---\left(15\right)$

然后将步骤(5)获得的几何图像矩向量Ψ_M代入式(13)，计算出未知三维Radon数据的投影 几何矩变换

进一步地，所述步骤(8)中投影几何矩逆变换的公式为：

${\widetilde{Rf}}_{{\bar{n}}^{\prime }}\left(\rho \right)={\Sigma }_{p=0}^M{\rho }^p{\tilde{L}}_p\left({\bar{n}}^{\prime }\right)---\left(17\right)$

其中，是未知三维Radon数据；

所述步骤(10)中三维Radon逆变换的公式为：

其中，是重建的CT图像。

采用上述技术方案带来的有益效果：

(1)本发明能够在减小扫描范围---即投影数据不完备条件下重建出符合临床诊断要 求、高质量的锥形束CT图像；

(2)由于本发明不涉及迭代运算，所以计算时间少于统计迭代方法；

(3)本发明可使用较少的锥形束投影数据重建出高质量的锥形束CT图像，故本发明 能在保证重建CT图像质量的前提下，有效地减少患者所受X射线辐射量。

附图说明

图1是本发明方法的基本流程示意图；

图2是锥形束投影的CT成像原理示意图；

图3是步骤(2)中单位向量的几何关系示意图；

图4是步骤(2)中直线t的几何关系示意图；

图5是使用传统的滤波反投影重建算法得到的重建图像；

图6是使用本发明得到的重建图像。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

基于几何图像矩的有限角度锥形束CT图像重建方法，如图1所示，包括如下步骤：

以三维Shepp-logan头模型为计算机仿真实验对象，尺寸为512*512*60，模拟CT机的X 射线源到旋转中心和探测器的距离分别为570mm和1040mm。锥形束CT设备旋转扫描范围 为[γ,360°-γ]，其中γ＝50°；角度采样值为1160；每个采样角对应672个探测器单元，探 测器单元的大小为1.407mm。

步骤(1)，获取有限角度扫描条件下的锥形束投影数据如附图2、3、4 所示，ρ为坐标原点O与特征点C之间的距离，直线t与线段OC_D垂直，且坐标原点O到直线t的 垂直距离为s；为从坐标原点O指向特征点C的单位向量， 由于CT设备未按常规沿轨道在[0°,360°)范围内旋转扫描一圈，而只是在有限角度 [γ,360°-γ]范围内旋转扫描。因此，对于已获得的锥形束投影数据其对应 的单位向量的参数值得注意的是，和角度范围 内的锥形束投影数据未知。

步骤(2)，利用Grangeat公式将已知锥形束投影数据转化为已知三维Radon数据关于ρ 的一阶导数。Grangeat公式为：

$R^{\prime }f\left(\rho \bar{n}\right)=\frac{1}{{\cos }^2\beta }\left(cos\alpha {\int }_{-\infty }^{\infty }\frac{\partial }{\partial p}{X}_{w}f\right(s\left(\rho \bar{n}\right),t)dt+\sin \alpha {\int }_{-\infty }^{\infty }\frac{\partial }{\partial q}{X}_{w}f(s\left(\rho \bar{n}\right),t\left)dt\right)---\left(1\right)$

其中，是已知三维Radon数据关于ρ的一阶导数；如附图2、3、4所示，cosβ＝SO/SC_D， SO是放射源S与坐标原点O之间的距离，C_D是放射线与探测器平面的交点，则SC_D是放射源 S与交点C_D之间的距离；ρ为坐标原点O与特征点C之间的距离；为从坐标原点O指向特征 点C的单位向量；是对步骤(1)中获取的不完备锥形束投影数据进 行加权计算后的投影数据值；所用的权值可优选为SO/SA，SA是放射源S与直线t上点的距 离；直线t与线段OC_D垂直，且坐标原点O到直线t的垂直距离为s，p,q分别为探测器平面的 横坐标轴和纵坐标轴，α为线段OC_D与探测器平面横坐标轴p的夹角；

然后对已知三维Radon数据关于ρ的一阶导数进行积分即可获得已知三维 Radon数据

步骤(3)，计算已知三维Radon数据的投影几何矩变换公式为：

$L_p\left(\bar{n}\right)={\int }_{-\infty }^{\infty }\rho ·{\mathrm{Rf}}_{\bar{n}}\left(\rho \right)d\rho ---\left(2\right)$

其中，为在方向上的P阶投影几何矩变换。

通过式(2)计算阶数P分别为[0,1,2,…,M]的一系列投影几何矩变换，M为所使用的几 何矩的最大阶数。在本实例中M＝20。然后，将这些投影几何矩变换值按一定顺序组成向 量

${\Phi }_M\left(\bar{n}\right)={[L_0\left(\bar{n}\right),L_1\left(\bar{n}\right),...L_M\left(\bar{n}\right)]}^T---\left(3\right)$

步骤(4)，针对已扫描角度范围，计算系数矩阵

其中，

步骤(5)，已知三维Radon数据的投影几何矩变换与几何图像矩之间关系的矩阵形式为：

${\Phi }_M\left(\bar{n}\right)={\Omega }_M\left(\bar{n}\right){\Psi }_M---\left(5\right)$

则利用公式(5)可得几何图像矩向量为：

${\Psi }_M={\Phi }_M\left(\bar{n}\right)/{\Omega }_M\left(\bar{n}\right)---\left(6\right)$

几何图像矩向量Ψ_M由一系列阶数不超过M的几何图像矩G_mnq组成，如下所示：

${\Psi }_M={[G^{{\left(0\right)}^T},G^{{\left(1\right)}^T},...G^{{\left(M\right)}^T}]}^T$

$G^{\left(k\right)}={[G_{k,0,M-k},G_{k-1,1,M-k},...,G_{1,k-1,M-k},G_{0,k,M-k}]}^T---\left(7\right)$

步骤(6)，针对未扫描角度范围，计算系数矩阵

其中，为未扫描角度范围内的向量，即或

步骤(7)，未知Radon数据的投影几何矩变换与几何图像矩之间关系的矩阵形式为：

${\tilde{\Phi }}_M\left({\bar{n}}^{\prime }\right)={\Omega }_M\left({\bar{n}}^{\prime }\right){\Psi }_M---\left(9\right)$

因为步骤(5)计算出的几何图像矩Ψ_M与均无关，所以可将步骤(5)获得的几何图 像矩Ψ_M以及步骤(6)获得的系数矩阵代入公式(9)，即可计算出未知Radon数据 的投影几何矩变换是一系列阶数不大于M的投影几何矩组成的向量 (0≤p≤M)。

步骤(8)，投影几何矩逆变换公式如下：

${\widetilde{Rf}}_{{\bar{n}}^{\prime }}\left(\rho \right)={\Sigma }_{p=0}^M{\rho }^p{\tilde{L}}_p\left({\bar{n}}^{\prime }\right)---\left(10\right)$

则利用公式(10)和步骤(7)获得的即可求得未知Radon数据

步骤(9)，将估计出的未知Radon数据按照扫描角度顺序和已知Radon数据合并，即可 获得补全的Radon数据。

步骤(10)，三维Radon逆变换公式如下：

利用三维Radon逆变换即可从补全的Radon数据中重建出CT图像

为了对比本发明所示方法的效果，使用传统的滤波反投影重建算法对不完备锥形束投 影数据进行重建，得到的重建图像如图5所示；而使用本发明方法获得的重建图像如图6 所示。比较图5和图6的重建结果，可以看出本发明方法所得结果中伪影相对较少，而且 头模型的轮廓更加完整和清晰。因此可以认为本发明方法能够在部分投影数据缺失的条件 重建出伪影更少、图像质量更好的结果。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照 本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明的保护范围之 内。