本发明涉及医学影像，特别涉及一种基于亚像素各向异性扩散的低剂量ct重建方法。
背景技术：
1、x射线ct因其在密度分辨力的卓越表现，已经广泛应用于临床医学影像诊断。然而ct扫描中过高的x射线辐射剂量会对人体造成辐射损伤、诱发恶性肿瘤等不可预测的伤害，所以在保证重建图像质量的前提下，如何最大限度地降低x射线剂量已成为医学ct成像领域的关键问题。
2、目前降低ct扫描过程中的管电流和x射线球管曝光时间是实现低剂量ct扫描最简单且最有效的方法。但是由于降低管电流和扫描时间，使得投影数据中含有大量的噪声，导致重建图像质量严重退化，难以满足临床诊断需要。
3、为了在保证图像质量的前提下大幅降低x射线辐射剂量，诸多低剂量ct图像重建方法相继提出，整体可分为三种策略进行：策略一为数据驱动方法，策略二为基于统计模型的迭代重建方法，策略三为基于低剂量ct投影数据恢复的图像重建方法。这三种策略均存在优点和缺点。
4、策略一是基于深度学习算法实现图像重建，通过学习大量的数据特征，不断优化模型，实现高精度的图像恢复。数据驱动技术因其高效的特征提取能力及强大的泛化能力，当前已在低剂量ct图像恢复中得到广泛研究。然而深度学习方法需要大量的数据进行训练，训练模型的复杂度较高，对数据的质量和规模都有很强的依赖性，且对计算资源的需求也会更高。同时，深度学习方法对于模型内部的运行机制缺乏可解释性，在医学图像处理中存在一定的限制。
5、策略二需要通过ct系统建模构建图像重建模型，优化目标函数实现图像重建。迭代重建算法充分考虑了投影数据的统计特性和ct系统的物理模型，而且易于加入先验信息约束，可以实现低剂量ct图像优质重建。然而重建同样大小的ct图像时，基于统计模型的迭代重建方法所花费的时间远远大于传统的解析重建方法。其原因在于基于统计模型的迭代重建方法，需要反复进行投影与反投影运算，使得迭代重建算法运算量大，图像重建时间大幅增加。
6、策略三首先对ct投影数据进行滤波，再采用经典的滤波反投影(filtered back-projection,fbp)算法对滤波后的投影数据进行重建。由于考虑到投影数据的噪声统计特性，这类方法具有高效且重建图像均匀性好等优点。然而由于传统的基于投影数据恢复的图像重建方法，在投影数据滤波过程中不可避免地会导致图像原有细节信息的丢失，从而导致相应ct图像分辨率的下降。
7、因此，针对现有技术不足，提供一种基于亚像素各向异性扩散的低剂量ct重建方法以解决现有技术不足甚为必要。
技术实现思路
1、本发明的目的在于避免现有技术的不足之处而提供一种基于亚像素各向异性扩散的低剂量ct重建方法。该基于亚像素各向异性扩散的低剂量ct重建方法解决了一般性单纯滤波降噪强度过大引起的图像细节丢失，在去除噪声和伪影的同时可以保持低剂量ct重建图像的边缘和结构细节特征。
2、本发明的上述目的通过以下技术措施实现：
3、提供一种基于亚像素各向异性扩散的低剂量ct重建方法，由如下步骤进行：
4、s1、获取投影数据，并对投影数据逐个数据点上的方差估计；
5、s2、构建投影数据的亚像素单元，然后计算亚像素单元的差分后，再通过求解各项异性扩散方程对投影数据去噪，得到去噪后投影数据，将根据s1得到的估方差对投影数据和去噪后投影数据进行加权平均，得到更新投影数据，然后对更新投影数据通过迭代最终得到待重建的投影数据；
6、s3、对s2得到的待重建的投影数据进行图像重建得到最终的重建ct图像。
7、优选的，上述s1具体为获取ct设备的系统参数和低剂量扫描协议下的投影数据yraw，并对投影数据yraw逐个数据点上的方差估计。
8、优选的，上述s2具体由如下步骤得到：
9、s2.1、令k＝0，将投影数据yraw定义为投影数据yk，进入s2.2；
10、s2.2、构建投影数据yk的亚像素单元然后计算亚像素单元差分，再求解各项异性扩散方程，对投影数据yk进行去噪，得到去噪后投影数据yfiltered，将s1得到的方差投影数据yraw及去噪后投影数据yfiltered进行加权平均得到更新投影数据yk+1，进入s2.3；
11、s2.3、判断是否到达第一迭代终止条件，当否时进入s2.4，当是时进入s2.5；
12、s2.4、令k+1＝k，使yk+1＝yk，返回s2.2；
13、s2.5、以当前的更新投影数据yk+1作为待重建投影数据。
14、所述第一迭代终止条件为其中||·||为欧几里得范数；或者所述第一迭代终止条件为k≥50。
15、优选的，上述s2.2具体由如下步骤得到：
16、s2.2.1、设定投影数据yk的亚像素单元距离hs，每两个相邻像素得到一个亚像素单元
17、s2.2.2、对于非边界像素的亚像素单元的值通过式(1)至式(4)的线性插值方法得到；
18、
19、
20、
21、
22、其中，和分别为的相邻像素，为投影数据yk在(i,j)处的值，i和j均为投影数据yk的像素索引，且i＝0,1,2...m，j＝0,1,2...n；和为周围的亚像素单元；
23、对于边界像素的亚像素单元的值，由投影数据对称计算获得，即以及
24、s2.2.3、将亚像素单元的差分通过一个离散的拉普拉斯算子表示，由式(5)至式(7)得到；
25、
26、
27、
28、其中，为在x方向的二阶差分，为在y方向的二阶差分，为的差分；
29、s2.2.4、将亚像素单元差分代替能量泛函的梯度，并将亚像素单元差分引入到各项异性扩散方程中，各向异性扩散方程在不同方向上的扩散过程由能量泛函表示，且所述能量泛函由式(8)表示；
30、
31、其中，ω为当前图像空间；
32、s2.2.5、所述能量泛函的最优解由式(9)的欧拉方程求解，得到去噪后投影数据yfiltered，进入s2.2.6；
33、
34、其中，t为当前时间步，div为散度，c(·)为扩散系数，且c(·)由式(10)表示；
35、
36、其中，f(·)为控制扩散速率的扩散系数方程，且f(·)由式(11)表示：
37、f(s)＝exp(-s2/ε2)……式(11)；
38、其中，ε为保护图像边缘的超参数，且ε＝3.0×10-3；
39、s2.2.6、将s1得到的方差投影数据yraw及s2.2.5得到的去噪后投影数据yfiltered，使用式(12)的加权平均滤波器进行加权平均，得到更新投影数据yk+1，进入s2.3；
40、
41、其中，β1为超参数。
42、本发明在所述s2.2.5中通过如下步骤对式(9)进行求解得到去噪后投影数据yfiltered，具体为：
43、s2.2.5.1、令t＝0，将投影数据yk定义为投影数据进入s2.2.5.2；
44、s2.2.5.2、将投影数据输入公式(13)得到更新投影数据进入s2.2.5.3；
45、
46、其中，λ为时间步长；
47、s2.2.5.3、判断是否到达第二迭代终止条件，当否时进入s2.2.5.4，当是时进入s2.2.5.5；
48、s2.2.5.4、令t+1＝t，使返回s2.2.5.2；
49、s2.2.5.5、停止迭代，将更新投影数据作为去噪后投影数据yfiltered，进入s2.2.6。
50、在所述s3中的图像重建为滤波反投影方法或者迭代重建方法。
51、优选的，上述系统参数为x射线入射光子强度i0或者系统电子噪声的方差中至少一种。
52、优选的，上述s1中的方差估计为基于小领域图像的局部方差估计或者基于投影数据噪声统计特性的方差估计。
53、第二迭代终止条件迭代次数为10次。
54、优选的，上述β1为1.0×10-3。
55、在所述s1中获得的投影数据为对数变换后的投影数据。
56、本发明的一种基于亚像素各向异性扩散的低剂量ct重建方法，由如下步骤进行：s1、获取投影数据，并对投影数据逐个数据点上的方差估计；s2、构建投影数据的亚像素单元，然后计算亚像素单元的差分后，再通过求解各项异性扩散方程对投影数据去噪，得到去噪后投影数据，将根据s1得到的估方差对投影数据和去噪后投影数据进行加权平均，得到更新投影数据，然后对更新投影数据通过迭代最终得到待重建的投影数据；s3、对s2得到的待重建的投影数据进行图像重建得到最终的重建ct图像。本发明的低剂量ct重建方法，通过引入亚像素单元能够通过组合多个整数像素单元的信息来提高图像分辨率和细节表现，因此可以在大幅减少图像噪声和伪影的同时较好地保持原有图像的分辨率，最终实现低剂量ct图像的优质重建，从而本发明解决了一般性单纯滤波降噪强度过大引起的图像细节丢失，在去除噪声和伪影的同时可以保持低剂量ct重建图像的边缘和结构细节特征。
1.一种基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于，由如下步骤进行：
2.根据权利要求1所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：所述s1具体为获取ct设备的系统参数和低剂量扫描协议下的投影数据yraw，并对投影数据yraw逐个数据点上的方差估计。
3.根据权利要求2所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于，所述s2具体由如下步骤得到：
4.根据权利要求3所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：所述第一迭代终止条件为其中·为欧几里得范数；
5.根据权利要求3所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于，所述s2.2具体由如下步骤得到：
6.根据权利要求5所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：在所述s2.2.5中通过如下步骤对式(9)进行求解得到去噪后投影数据yfiltered，具体为：
7.根据权利要求1至6任意一项所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：在所述s3中的图像重建为滤波反投影方法或者迭代重建方法。
8.根据权利要求1至6任意一项所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：所述系统参数为x射线入射光子强度i0或者系统电子噪声的方差中至少一种。
9.根据权利要求1至6任意一项所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：所述s1中的方差估计为基于小领域图像的局部方差估计或者基于投影数据噪声统计特性的方差估计。
10.根据权利要求6所述的基于亚像素各向异性扩散的低剂量ct重建方法，其特征在于：所述第二迭代终止条件迭代次数为10次；
