基于噪声生成机制与数据驱动紧框架的CT投影去噪重建方法及装置与流程

本发明属于ct图像重建技术领域，特别涉及一种基于噪声生成机制与数据驱动紧框架的ct投影去噪重建方法及装置。

背景技术：

计算机断层成像(computedtomography,ct)利用x射线穿透物体的透射信息，通过重建算法计算物体内部衰减特性分布。作为20世纪后期最伟大的科技成果之一，ct技术可以通过非接触、非破坏的方式高效率地获取物体内部结构信息。1972年，第一台ct设备由工程师hounsfield成功研发，随后在医学领域获取了快速发展，掀起了临床疾病诊断的革命。ct可以提供病人清晰准确的结构图像，有利于医师高效准确地辨别病变，目前已成为临床疾病诊断不可或缺的工具之一。

随着ct技术在临床上的广泛应用，过多辐射剂量所造成的危害逐渐成为人们关注的焦点。当辐射剂量超过安全范围时，会影响人体的身体健康，增加癌症、白血病等疾病的患病风险。1990年，naidich等人首次提出低剂量ct(low-dosecomputedtomography,ldct)的概念，引起了ct领域研究者的广泛关注。目前，ldct的实现主要有两种方式，一是通过降低x射线源的管电流来降低单角度投影采集时的x射线强度；二是减少投影数量，即稀疏角度问题。由于实际临床扫描速度较快，稀疏采样不易实现，而降低管电流可操作性强，简便有效，目前已成为临床上实现ldct的主要手段。但降低管电流时，由于探测器接收到的光子数大量减少，采集得到的投影往往被光子统计噪声和探测系统热噪声严重污染，是ct图像质量严重下降的本质原因。因此，亟需一种保证去噪以及细节保留能力的低剂量ct投影域重建方法，以解决ct临床中过多辐射剂量的危害。

技术实现要素：

为此，本发明提供一种基于噪声生成机制与数据驱动紧框架的ct投影去噪重建方法及装置，克服低剂量ct中噪声的严重影响，提高ct图像重建质量。

按照本发明所提供的设计方案，一种基于噪声生成机制与数据驱动紧框架的ct投影去噪重建方法，包含如下内容：

刻画投影数据噪声统计特性，并挖掘投影数据的稀疏先验信息；

结合噪声统计特性和稀疏先验信息，建立最大后验概率模型；

对最大后验概率模型进行迭代求解，获取去噪后的正弦图数据；

对去噪后的正弦图数据进行重建，得到重建后的ct图像。

上述的，刻画投影数据噪声统计特性中，通过分析投影数据噪声本质来源，考虑光子统计噪声和系统热噪声统计特性，建立噪声生成机制，通过该噪声生成机制对投影数据的统计特性进行刻画。

优选的，挖掘投影数据的稀疏先验信息中，将稀疏先验信息嵌入噪声生成机制，并通过数据驱动紧框架挖掘投影数据的稀疏先验信息。

上述的，建立最大后验概率模型，包含如下内容：首先，采用基于小波变换的数据驱动紧框架对投影数据进行稀疏编码，获取投影数据和先验信息的联合分布；然后，结合噪声生成机制，获取后验概率分布；基于该后延概率分布，采用map框架估计正弦图，并得到最大后验概率模型。

上述的，对最大后验概率模型进行迭代求解中，将模型求解等价为带约束的优化问题；并利用增广拉格朗日函数将带约束的优化问题转化为无约束的优化问题；采用交替方向承子admm法将该无约束的优化问题进行问题分解；根据预先设定的循环迭代条件，对分解后的问题进行迭代求解，得到去噪后的正弦图数据。

优选的，用交替方向承子admm法将该无约束的优化问题进行问题分解中，依据接收射线量子、正弦图及尺度参数进行问题分解，获取若干子问题；然后，根据预先设定的循环迭代条件，对分解后的该若干子问题进行迭代求解，得到去噪后的正弦图数据。

优选的，对分解后的若干子问题迭代求解过程，包含如下内容：初始化求解参数，获取投影数据和正弦图；依次对子问题、求解参数和小波基函数进行参数更新求解，并判断是否达到预设最大迭代次数，若是，则判定当前更新求解过程的正弦图数据即为去噪后的正弦图数据，否则，返回参数更新求解步骤继续执行。

上述的，采用经典解析fbp算法对去噪后的正弦图数据进行重建，得到重建后的ct图像。

一种基于噪声生成机制与数据驱动紧框架的ct投影去噪重建装置，包含：数据处理模块、模型构建模块、模型求解模块和图像重建模块，其中，

数据处理模块，用于刻画投影数据噪声统计特性，并挖掘投影数据的稀疏先验信息；

模型构建模块，用于结合噪声统计特性和稀疏先验信息，建立最大后验概率模型；

模型求解模块，用于对最大后验概率模型进行迭代求解，获取去噪后的正弦图数据；

图像重建模块，用于对去噪后的正弦图数据进行重建，得到重建后的ct图像。

上述的装置中，模型构建模块包含：稀疏编码单元、概率分布计算单元和模型获取单元，其中，

稀疏编码单元，用于采用基于小波变换的数据驱动紧框架对投影数据进行稀疏编码，获取投影数据和先验信息的联合分布；

概率分布计算单元，用于结合噪声生成机制，获取后验概率分布；

模型获取单元，用于基于后延概率分布，采用map框架估计正弦图，并得到最大后验概率模型。

本发明的有益效果：

本发明中，首先通过刻画投影数据噪声统计特性，并挖掘投影数据的稀疏先验信息；结合两者数据建立最大后验概率模型；然后对模型迭代求解，获取去噪后的正弦图数据；通过重建获取ct图像；克服了传统低剂量ct重建算法因采用近似的方式刻画投影数据的噪声导致重建图像易残留严重的噪声和伪影的问题；并进一步引入噪声生成机制来充分刻画投影数据固有的“泊松+高斯”统计特性，同时利用基于学习策略的数据驱动紧框架对投影数据进行稀疏信息挖掘，利用admm算法对模型进行高效求解，能够在尽可能保留投影数据有用信息的同时去除噪声，进而有效抑制了重建图像的噪声和伪影，在提升低剂量ct重建图像质量方面具有较高的应用价值。

附图说明：

图1为实施例中方法流程示意图；

图2为实施例中计算过程示意图；

图3为实施例中装置示意图；

图4为实施例中模型构建模块示意图；

图5为实施例中仿真实验重建图像对比图。

具体实施方式：

为使本发明的目的、技术方案和优点更加清楚、明白，下面结合附图和技术方案对本发明作进一步详细的说明。

基于投影预处理的ldct重建算法的基本思路是首先对ldct系统直接采集的投影数据进行预处理以抑制投影域噪声，而后采用fbp算法对去噪后的投影数据进行重建，从而改善ldct图像的质量。现阶段国内外学者对该类算法的研究主要分为两种方式：一种是通过改进滤波器的性能以实现去噪效果的改善；另一种是基于投影数据的统计特性建立统计迭代模型对投影数据进行噪声抑制。该种方法将投影数据去噪过程看作优化函数问题，在投影数据统计特性的基础上，建立统计迭代模型进行去噪，由于这种类方法更具有针对性，因此受到了国内外学者的广泛关注并在实际应用中取得了较好的效果。通过对ldct多次扫描所得投影数据的统计特性进行探究，发现对采集的投影数据进行校正和对数变换后，其噪声服从一种非平稳高斯分布，且方差与均值满足非线性关系。为此，本发明实施例提出基于迁移学习的fmri视觉编码模型。所谓迁移学习，通俗来说就是将从源领域学习到的东西应用到新的领域。迁移学习可以解决小数据问题，能够将用于大数据的模型迁移到小数据上。本发明实施例中，参见图1所示，提供一种基于噪声生成机制与数据驱动紧框架的ct投影去噪重建方法，包含如下内容：

s101、刻画投影数据噪声统计特性，并挖掘投影数据的稀疏先验信息；

s102、结合噪声统计特性和稀疏先验信息，建立最大后验概率模型；

s103、对最大后验概率模型进行迭代求解，获取去噪后的正弦图数据；

s104、对去噪后的正弦图数据进行重建，得到重建后的ct图像。

首先通过刻画投影数据噪声统计特性，并挖掘投影数据的稀疏先验信息；结合两者数据建立最大后验概率模型；然后对模型迭代求解，获取去噪后的正弦图数据；通过重建获取ct图像；克服了传统低剂量ct重建算法因采用近似的方式刻画投影数据的噪声导致重建图像易残留严重的噪声和伪影的问题，提升图像重建质量。

本发明实施例中，刻画投影数据噪声统计特性中，通过分析投影数据噪声本质来源，考虑光子统计噪声和系统热噪声统计特性，建立噪声生成机制，通过该噪声生成机制对投影数据的统计特性进行刻画。优选的，挖掘投影数据的稀疏先验信息中，将稀疏先验信息嵌入噪声生成机制，并通过数据驱动紧框架挖掘投影数据的稀疏先验信息。通过引入噪声生成机制以充分刻画投影数据的统计特性并设计正则化项充分挖掘先验信息进而提高模型的去噪以及细节保留能力。

本发明另一个实施例，建立最大后验概率模型中，采用基于小波变换的数据驱动紧框架对投影数据进行稀疏编码，获取投影数据和先验信息的联合分布；然后，结合噪声生成机制，获取后验概率分布；基于该后延概率分布，采用map框架估计正弦图，并得到最大后验概率模型。

上述的，对最大后验概率模型进行迭代求解中，将模型求解等价为带约束的优化问题；并利用增广拉格朗日函数将带约束的优化问题转化为无约束的优化问题；采用交替方向承子admm法将该无约束的优化问题进行问题分解；根据预先设定的循环迭代条件，对分解后的问题进行迭代求解，得到去噪后的正弦图数据。用交替方向承子admm法将该无约束的优化问题进行问题分解中，依据接收射线量子、正弦图及尺度参数进行问题分解，获取若干子问题；然后，根据预先设定的循环迭代条件，对分解后的该若干子问题进行迭代求解，得到去噪后的正弦图数据。优选的，对分解后的若干子问题迭代求解过程，包含如下内容：初始化求解参数，获取投影数据和正弦图；依次对子问题、求解参数和小波基函数进行参数更新求解，并判断是否达到预设最大迭代次数，若是，则判定当前更新求解过程的正弦图数据即为去噪后的正弦图数据，否则，返回参数更新求解步骤继续执行。

投影重建技术中，fbp算法是在傅立叶变换理论基础之上的一种空域处理技术。它的特点是在反投影前将每一个采集投影角度下的投影进行卷积处理，从而改善点扩散函数引起的形状伪影，重建的图像质量较好。本发明另一个实施例中，采用经典解析滤波反投影fbp算法对去噪后的正弦图数据进行重建，得到重建后的ct图像。

本发明实施例中具体计算过程可设计为如下内容：

参见图2所示，步骤1中：建立噪声生成机制，充分刻画投影数据的噪声统计特性：

在理想情况下，由beer-lambert定律可知，投影测量数据的表达式为：

其中为一维向量，表示正弦图，即测量数据对数变换后的数据，n表示投影路径数目。i0i(i＝1,…n)表示沿投影路径i的入射x射线的强度。

由上述噪声来源分析可知，测量数据f不可避免的存在噪声，可采用如下模型进行描述：

f＝q+ε,(2)

其中表示探测器接收到的射线的量子，表示系统热噪声，分别表示实数集合、自然数集合，一般假设该噪声服从高斯分布，可用下式表示：

其中σ²为噪声的方差。

由于x射线在发射过程中会不可避免地出现量子波动，因此探测器接收到的量子可看作理想射线被扰动的结果，泊松分布可以较好刻画该过程：

其中o为理想情况下探测器接收到的射线量子，即射线量子的平均数量。令y为理想的正弦图，由beer-lambert定律可得：

结合式(4)和(5)可得：

由式(3)可知：

结合式(6)和(7)可得：

通过分析投影数据噪声的本质来源，式(8)综合考虑了光子统计噪声和系统热噪声的统计特性，充分利用了投影数据的统计特性，相比于传统刻画方式，能更充分地描述隐含在正弦图中的噪声生成机制。将在该噪声机制框架的基础上，设计针对性去噪算法对低剂量投影数据进行去噪。

步骤2中：采用hessian惩罚项将稀疏先验模型嵌入到噪声生成机制模型中。并选用基于学习策略的数据驱动紧框架来挖掘投影数据的稀疏先验信息，建立基于噪声生成机制与数据驱动紧框架的最大后验概率模型：

对于稀疏先验信息的引入，本文采用hessian惩罚项将稀疏先验模型嵌入到统计模型中。通常来讲，可以采用如下公式对数据进行稀疏编码：

其中，为指数幂分布(也称为超拉普拉斯分布)，b为尺度参数，f:为线性算子，cf为与算子f相关的常数。通过采用基于小波变换的数据驱动紧框架对投影数据进行稀疏编码。令w表示小波算子，则f(y)＝wy。假设参数b的先验分布为p(b)(通常为常数函数)，则由式(9)可得数据y和b的联合分布p(y,b)：

结合式(8)和式(10)可得到如下后验概率分布：

其中p(f)为与观测值f相关的常数。基于该后验概率分布，我们可以采用map框架来估计正弦图y，即解决如下优化问题：

步骤3中：采用交替方向乘子法设计高效求解算法：令z＝wy，则式(12)可改写为：

s.t.wy＝z.

对于该模型，可等价为求解如下带约束的优化问题：

s.t.wy＝z.

利用增广lagrangian函数法将带约束的优化问题(14)转化为如下无约束的优化问题：

其中，为lagrangian乘子，μ为非负常数。

采用admm方法将lμ(q,y,z,b,λ)原问题分解为四个子问题，实现对上述优化模型的求解，分解得到的四个子问题分别为l(q)、l(y)、l(z)、l(b)。

(1)子问题l(q)可表示为：

对于每一个qi可采用下式进行求解：

(2)子问题l(y)可表示为：

对于该优化问题，通过采用加速一阶方法对其进行求解：给定采用二次逼近的方法对式(18)的第二项进行近似，可得：

其中⊙表示逐元素相乘，参数τ可设置为i0的最大值。对式(19)求导并令其导函数为0，则可得到如下公式：

由数据驱动紧框架的特性可得w^tw＝i，此时y的求解公式为：

(3)子问题l(z)可表示为：

该式的局部最优解的闭式解为：

其中sλ(·)为门限算子，定义如下：

其中

(4)子问题l(b)可表示为：

对上式求导并令其导函数为0，可得：

步骤4：实现并运行算法，获得较低噪声的重建图像：现给出基于噪声生成机制与数据驱动紧框架的投影去噪算法，获得低剂量ct重建图像的具体方法如下：

输入：i0i，初始μ,b,ρ,λ＝0和y＝max{ln(i0)-ln(f),0}，执行如下迭代过程：(1)更新q子问题，对于每一个qi，公式如下：

(2)更新y子问题，公式如下：

(3)更新z子问题，公式如下：

(4)更新b子问题，公式如下：

(5)参数更新，公式如下：

(6)更新小波基函数，公式如下：

重组操作：

{yy,zz,h}←{y^k+1,z^k+1,wl},(32)

更新小波基函数：

mσn^t＝zzyy^t,hl+1←mn^t,(33)

更新变量z：

zz←γλ((hl+1)^tyy),(34)

变量替换：

{z,w}←{zz,hl+1},(35)

(7)迭代轮数增加：k←k+1。若“未达到迭代设定的最大次数”则返回第(1)步，否则停止，迭代达到最大轮数退出时，y⁽ⁿ⁾即输出去噪后的正弦图数据，然后采用经典解析算法fbp对去噪后正弦图进行重建即可得到重建后的ct图像。

基于上述的重建方法，本发明实施例还提供一种基于噪声生成机制与数据驱动紧框架的ct投影去噪重建装置，参见图3所示，包含：数据处理模块101、模型构建模块102、模型求解模块103和图像重建模块104，其中，

数据处理模块101，用于刻画投影数据噪声统计特性，并挖掘投影数据的稀疏先验信息；

模型构建模块102，用于结合噪声统计特性和稀疏先验信息，建立最大后验概率模型；

模型求解模块103，用于对最大后验概率模型进行迭代求解，获取去噪后的正弦图数据；

图像重建模块104，用于对去噪后的正弦图数据进行重建，得到重建后的ct图像。

上述的装置中，模型构建模块102包含：稀疏编码单元2001、概率分布计算单元2002和模型获取单元2003，其中，

稀疏编码单元2001，用于采用基于小波变换的数据驱动紧框架对投影数据进行稀疏编码，获取投影数据和先验信息的联合分布；

概率分布计算单元2002，用于结合噪声生成机制，获取后验概率分布；

模型获取单元2003，用于基于后延概率分布，采用map框架估计正弦图，并得到最大后验概率模型。

为进一步说明本发明的有效性，本发明实施例中还通过仿真实验对本发明的有效性做进一步验证，参见图5所示，(a)为参考图像，(b)为低剂量ct重建图像，(c)为本发明实施例技术方案的重建图像；通过图示也可看出，本发明中方案能够在尽可能不破坏投影数据有用信息的基础上有效抑制噪声，对重建图像的噪声和伪影具有较好的去除效果，在提升低剂量ct重建图像质量方面具有较高的应用价值

除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对步骤、数字表达式和数值并不限制本发明的范围。

基于上述的方法，本发明实施例还提供一种服务器，包括：一个或多个处理器；存储装置，用于存储一个或多个程序，当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现上述的方法。

基于上述的方法，本发明实施例还提供一种计算机可读介质，其上存储有计算机程序，其中，该程序被处理器执行时实现上述的方法。

本发明实施例所提供的装置，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述方法实施例中相应内容。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

在这里示出和描述的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制，因此，示例性实施例的其他示例可以具有不同的值。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read-onlymemory)、随机存取存储器(ram，randomaccessmemory)、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。