技术领域
[0001] 本
发明涉及一种高能X射线图像非盲去模糊方法，属于
图像处理技术领域。
背景技术
[0002] 高能X射线成像技术是研究客体内部结构的一种重要手段，是核武器研究的主要工具，具有成像直观，灵敏度高，对物体不破坏、无污染，适用材料广泛等优点，在医学、工业、材料和国防军事等很多领域中都有极为广泛的应用。由于高能X射线成像系统本身很复杂，成像过程会受到模糊、噪声等影响，致使图像
质量较差。为了从模糊的高能X射线图像中获得更佳的诊断结果，发展了各种图像处理方法。传统的解析方法、
迭代方法以及统计方法等在医学影像诊断、工业无损探测等领域得到了广泛的应用，但却难以满足高能X射线成像这类复杂场合下的高维反演要求。
[0003] 针对高能X射线图像，求解的问题属于高维反演，面对的主要难题是消除模糊及噪声的影响。目前，主要反演方法有
滤波反投影法(FBP)、约束共轭梯度法(CCG)及基于不同正则化技术的反演方法等。上述方法中FBP反演图像
分辨率低，对噪声的抑制能
力较差；CCG主要考虑了非负约束和光滑性约束条件，还没有对系统模糊影响边界恢复质量的问题研究针对性的约束条件；基于正则化技术的反演方法主要将正则项作为惩罚项，通过正则化参数来调节数据保真项与正则项之间的矛盾，没有考虑引入先验约束条件。
发明内容
[0004] 本发明所要解决的技术问题是克服
现有技术的高能X射线图像去模糊的
进程过慢，高能X射线图像的清晰度不够的
缺陷，提供一种高能X射线图像非盲去模糊方法。
[0005] 为解决上述技术问题，本发明提供一种高能X射线图像非盲去模糊方法，获取高能X射线模糊图像y和预先已给定的模糊核K；
[0006] 将高能X射线模糊图像y和预先给定的模糊核K作为输入至预先构建的频域上的最大后验概率模型，并设定该模型的
采样次数Nsamps和老化的样本数Nb，初始化λ0和δ0，得到变量x，λ和δ的后验
密度函数；
[0007] 运用Gibbs抽样从变量x，λ和δ的后验密度函数中提取样本，动态构造
马尔可夫链，当各变量的各条马尔可夫链达到平稳后，进行抽样，得到满足同分布的各变量，计算平稳后的各变量样本均值，得到初步去模糊图像xa、正则化参数τ和噪声标准差1/λ的估计量，其中变量x表示待估计的清晰图像，λ表示噪声
精度参数，δ表示先验精度参数；
[0008] 利用迭代引导滤波对初步去模糊图像xa进行噪声的去除，并通过噪声标准差1/λ进行迭代次数的判断，得到滤波图像xfilter；
[0009] 基于x的统计样本的贝叶斯置信区间对初步去模糊图像xa的负值
像素进行非负插值处理，并在滤波图像xfilter上，用处理后的值替换掉对应于xa负值
位置上的值，同时抑制滤波图像xfilter的
振铃效应，得到平滑图像
[0010] 基于稀疏引导正则项约束的子目标函数，将处理后的平滑图像 以及估计的正则化参数τ作为该目标函数的暖启动，获得最终平滑且边缘信息丰富的高能X射线去模糊图像。
[0011] 进一步的，所述高能X射线模糊图像y利用第一类弗雷德霍姆积分方程来表示，经过离散化后得到一个线性方程，其表达式为：
[0012] y＝Kx+η(1)
[0013] 其中，x∈RN为未知的清晰图像；y∈RM为高能X射线模糊图像；K∈RM×N为模糊核数值离散后得到的正向矩阵；η∈RM为加性噪声，RN和RM×N分别表示N维和M×N维的实数空间，xN M M×N∈R，y∈R表示x,y分别为N×1和M×1的列向量，K∈R 表示为M×N的矩阵。
[0014] 进一步的，所述图像非盲去模糊是一个病态逆问题，常常采用正则化方法进行建模，设计如下的目标函数对去模糊问题进行求解：
[0015]
[0016] 其中，第一项为数据保真项；第二项为正则项；xτ为恢复的清晰图；τ是 范数和稀疏引导的 范数的统一正则化参数，用来控制保真项与正则项之间的加权比例；L为正则化矩阵；
[0017] 为了避免人为选择正则化参数，利用统计方法，将目标函数转化为以稀疏引导的范数作为正则项约束的目标函数，为：
[0018]
[0019] 其中，
[0020]
[0021] 式(3)为以图像 范数作为正则项约束的子目标函数。
[0022] 式(3)引入 范数使得图像更加平滑，同时抑制了图像中的噪声，得到经过式(3)处理的平滑图像后，式(4)引入稀疏引导的 范数能更好地保留图像的边缘信息。
[0023] 进一步的，所述最大后验概率模型表示为：
[0024]
[0025] 其中，N(,)表示正态分布，Γcov为x|y,λ,δ所服从的正态分布的方差，T为转置操作，Γprior为负拉普拉斯算子的离散化。
[0026] 进一步的，采用迭代引导滤波进行噪声的去除，并通过噪声标准差1/λ根据进行迭代次数的判断，最后得到滤波图像xfilter，表示为：
[0027]
[0028] 其中，i和j为像素标签，itr为迭代的次数，当itr＝1时，xaitr就是xa；当itr＞1时，上一次的滤波图像则作为下一次滤波的引导图像，Gij为滤波核函数，将Gij定义为：
[0029]
[0030] 其中，wk为第k个核函数窗口， 为窗口内的像素总数，μk和 是引导图像在窗口内的均值和方差，ε为平滑因子。
[0031] 进一步的，得到图像 的过程为：
[0032] 将初步去模糊图像xa的95％的置信区间设为 xa的每个像素都有置信区间；
[0033] 计算 其中，Nnv为负值的数量，(i,j)表示像素，Loc_nv表示图像xa中所有负值所在的像素位置的集合；
[0034] 利用噪声标准差1/λ根据 进行迭代次数的判断，并通过式(18)对 和 进行迭代引导滤波，u和l分别表示置信区间的上边界和下边界；
[0035] 针对Loc_nv的Nnv个负值进行非负插值，具体步骤如下：
[0036] ①令t＝1:Nnv，并定义ψt为图像xa中以 为中心的大小为r×r的窗口，其中，表示图像xa在Loc_nv(t)处的负值像素，再对ψt中的非负像素按升序排序，记为ψ't；
[0037] ②如果 的95％置信区间的上边界 小于0或者 大于0但小于ψ't窗口中的最小像素值时，计算 得到ψ't窗口中符合要求的像素个数 其
中， 表示正向取整操作；ρ的作用是控制在非负插值中起作用的像素的百分比，再计算表示ψ't窗口中的像素集合，用 代替xfilter中对应位置
上的值，得到图像
[0038] ③如果不满足步骤②中的条件，则计算ψ't窗口中的像素值小于等于下边界的像素个数Nu或者计算ψ't窗口中的像素值小于等于上边界 且大于等于下边界 的像素个数Nu，再计算 用 代替xfilter中对应
位置上的值，得到图像
[0039] 进一步的，将处理后的平滑图像 以及估计的正则化参数τ作为引入了稀疏引导的范数的子目标函数的暖启动，最终获得平滑且边缘细节丰富的高能X射线去模糊图像。
[0040] 进一步的，利用GPSR-BB
算法将包含 范数正则项的子目标函数转化为有界约束的二次规划方程进行求解。
[0041] 引入稀疏引导的 范数作为子目标函数的正则项，能在保持图像平滑的同时，更好地保留边缘信息。针对包含 范数正则项的子目标函数，GPSR-BB算法可将其转化为有界约束的二次规划方程进行求解。于是，利用GPSR-BB算法对 进行非盲去模糊，得到最终的高能X射线去模糊图像。
[0042] 本发明所达到的有益效果：
[0043] 本发明可以
加速高能X射线图像的去模糊进程，更好地去除图像中的噪声和振铃，在保持平滑度的同时获得更加精确的去模糊结果及其不确定度的统计量，而不会产生额外的噪声。
具体实施方式
[0044] 下面对本发明作进一步描述。以下
实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
[0045] 一种高能X射线图像非盲去模糊方法，具体步骤如下：
[0046] 步骤1：获取高能X射线模糊图像y，构造去模糊的目标函数；
[0047] 图像模糊过程可利用第一类弗雷德霍姆积分方程来表示，经过离散化后得到一个线性方程，其表达式为：
[0048] y＝Kx+η(1)
[0049] 其中，x∈RN为未知的清晰图像；y∈RM为观测的模糊图像；K∈RM×N为模糊核数值离散后得到的正向矩阵，又称为点扩散函数(PSF)；η∈RM为加性噪声，在本发明中，η是一个均值为0，方差为1/λ的高斯随机变量，RN和RM×N分别表示N维和M×N维的实数空间，x∈RN，y∈RMM×N说明x,y分别为N×1和M×1的列向量，K∈R 说明K为M×N的矩阵。
[0050] 图像非盲去模糊是假设模糊图像和模糊核都已给出，估计清晰图像。在图像的恢复处理中，图像非盲去模糊是一个病态逆问题，常常采用正则化方法进行建模。本发明设计如下的目标函数对去模糊问题进行求解：
[0051]
[0052] 其中，第一项为数据保真项；第二项为正则项；xτ为恢复的清晰图；τ是二次范数和稀疏引导范数的统一正则化参数，用来控制保真项与正则项之间的加权比例；L为正则化矩阵。
[0053] 为了避免人为选择正则化参数，利用统计方法，将目标函数转化为以稀疏引导的范数作为正则项约束的目标函数，为：
[0054]
[0055] 其中，
[0056]
[0057] 式(3)为以图像 范数作为正则项约束的子目标函数。
[0058] 式(3)引入 范数使得图像更加平滑，同时抑制了图像中的噪声。得到经过式(3)处理的平滑图像后，式(4)引入稀疏引导的 范数能更好地保留图像的边缘信息。
[0059] 步骤2：基于给定的模糊核K和高能X射线模糊图y，采用贝叶斯模型建立最大后验概率(MAP)模型，在频域上运用Gibbs抽样动态构造马尔可夫链，进行快速非盲去模糊，进而获得图像的初步去模糊结果xa，正则化参数τ以及噪声标准差1/λ的估计量；所述Gibbs抽样是目前MCMC方法中一种有效的高维抽样方法，其基本思想就是从各变量的条件概率分布中，对各变量进行抽样。其中就包含了对x的抽样，当各条马尔可夫链达到平稳后，再进行抽样得到的x是满足同分布的，计算平稳后的样本均值，可以认为其就是对原始数据的一个有效估计；
[0060] 最大后验概率模型是根据经验数据对难以观察的量进行估计。与最大似然模型类似，但最大的不同在于，最大后验概率模型融入了所要估计量的先验分布在其中，所以最大后验概率模型可看作是规则化的最大似然估计模型。应用贝叶斯理论建立最大后验概率模型，定义采样的后验密度函数，去实现x|y的概率最大化。首先，通过式(1)定义似然函数：
[0061]
[0062] 其中，T为转置操作，λ表示噪声精度参数，∝表示成比例。接着，把对未知x的不确定性，以及对x的某些性质的先验知识进行编码，并假设x～N(0,δ-1Γprior)，得到x|δ的先验概率密度函数：
[0063]
[0064] 其中，Γprior为负拉普拉斯算子的离散化，并且， δ为先验精度参数。因此，以y,λ,δ为条件的x的后验概率密度函数为：
[0065] p(x|y,λ,δ)∝p(y|x,λ)p(x|δ)(7)
[0066] 最大后验概率估计是使得p(x|y,λ,δ)最大化，也就是说最小化-lnp(x|y,λ,δ)。在式(3)中，δ/λ与正则化参数τ是等价的。为了可以不事先选择最优参数，本发明中将δ和λ当作未知的，对它们施加独立共轭的伽马分布，为它们设置初始值，取δ＝0.005，λ＝7，并定义δ和λ的联合概率密度函数：
[0067]
[0068] 其中，αλ、αδ分别表示关于λ和δ的伽马分布的形状参数，βλ、βδ分别表示关于λ和δ的伽马分布的逆尺度参数，p(λ)和p(δ)分别表示参数λ和δ的共轭超先验。本发明中将αλ，αδ设置为1，βλ，βδ设置为10-4，使得关于δ和λ的伽马分布的均值和方差比正常值高几个数量级。接着，运用Gibbs方法在采样过程中估计出δ和λ这两个精度参数，并定义联合后验概率密度函数为：
[0069]
[0070] 根据联合后验概率密度，推导出如下的全条件概率密度函数，并将MCMC采样应用于样本进行关于x,λ和δ的统计特性分析。最后，从样本的均值和方差中获得图像的初步去模糊结果xa以及正则化参数τ和噪声标准差1/λ的估计量。
[0071]
[0072]
[0073]
[0074] 在时域上，运用MCMC方法进行非盲去模糊时，像y∈Rm×n这样的图像需要经过y＝vec(y)以向量y∈RM＝m×n的形式参与计算，并且将模糊核k扩展成正向矩阵K∈RM×N，所以每一次采样都需要巨大的计算成本。在计算机
断层扫描中，经常在时域上进行计算。但在图像去模糊问题中，本发明在频域上加速了去模糊问题的求解，频域上的最大后验概率模型为：
[0075]
[0076] 其中，N(,)表示正态分布，Γcov为x|y,λ,δ所服从的正态分布的方差。
[0077] 式(13)的有效闭合解为：
[0078]
[0079]
[0080] 其中， 和 分别代表快速傅里叶变换(FFT)和反快速傅里叶变换(IFFT)；代表复共轭算子；γ1和γ2为正态分布的随机矩阵； 的定义是建立在高斯马尔可夫随机场上的，即假设，全条件概率密度函数均为高斯分布，与 的联合概率密度函数也为高斯分布，则称 为高斯马尔可夫随机场。在周期边界的条件下，与 相乘相当于与一个大小为m×n的矩阵l卷积。
[0081]
[0082] 基于式(13)，即频域上的最大后验概率模型，将高能X射线模糊图y和预先给定的模糊核K作为输入，并设定采样次数Nsamps和老化的样本数Nb，本发明中，取Nsamps＝200，Nb为Nsamps的10％，然后初始化λ0和δ0。从x，λ和δ的后验密度函数中提取样本，即对各变量进行Gibbs抽样，动态构造马尔可夫链。当各变量的各条马尔可夫链达到平稳后，再进行抽样得到的各变量是满足同分布的。计算平稳后的样本均值，可认为其就是对原始数据的一个有效估计，于是可得到初步去模糊图像xa以及正则化参数τ和噪声标准差1/λ的估计量。具体抽样步骤如下：
[0083] ①通过式(14)计算xt；
[0084] ②计算
[0085] ③计算
[0086] ④令t＝t+1并返回步骤①，直到达到采样次数Nsamps；
[0087] ⑤从抽样样本中估计出x,λ和δ，并根据 计算出初步去模糊图像xa。
[0088] 步骤3：利用迭代引导滤波对步骤2得到的初步去模糊图像xa进行噪声的去除，并通过步骤2估计的噪声标准差1/λ进行迭代次数的判断，得到滤波图像xfilter；
[0089] 样本重建后的初步去模糊图像xa仍然存在噪声，为了精细图像，采用迭代引导滤波来实现噪声的去除，并通过步骤2估计的噪声标准差1/λ根据 来进行迭代次数的判断，最后得到滤波图像xfilter。
[0090] 图像引导滤波是一个线性移可变的滤波过程，包括引导图像，输入图像和输出图像。其中引导图像是需要根据具体应用事先设定的，也可以直接取为输入图像。本发明中，迭代引导滤波的输入图像和引导图像均为xa，输出图像为xfilter。对于xfilter的第i个像素而言，其计算方法可表示为：
[0091]
[0092] 其中，i和j为像素标签，itr为迭代的次数，当itr＝1时，xaitr就是xa；当itr＞1时，上一次的滤波图像则作为下一次滤波的引导图像。Gij为滤波核函数，将Gij定义为：
[0093]
[0094] 其中，wk为第k个核函数窗口； 为窗口内的像素总数；μk和 是引导图像在窗口内的均值和方差，ε为平滑因子。
[0095] 步骤4：基于MCMC统计样本95％的置信区间，对步骤2得到的初步去模糊图像xa的负值像素进行非负处理，并在图像xfilter上，用处理后的值替换掉对应于xa负值位置上的值，同时抑制图像xfilter的振铃效应，得到图像 本发明中，95％的置信图像由每个像素点上样本的0.025和0.975经验分位数计算得到。
[0096] 图像的像素均不小于0。但步骤2的Gibbs抽样会产生负值，由于负值的影响，图像xfilter会存在振铃效应。为了抑制图像中的振铃，本发明利用非负插值来剔除负值。
[0097] 首先，将步骤2得到初步去模糊图像xa的95％的置信区间设为 u和l分别表示置信区间的上边界和下边界，xa的每个像素都有置信区间；
[0098] 接着，计算 其中，Nnv为负值的数量，(i,j)表示像素，Loc_nv表示图像xa中所有负值所在的像素位置的集合；
[0099] 然后，利用步骤2估计的噪声标准差1/λ根据 进行迭代次数的判断，并通过(18)对 和 进行迭代引导滤波；
[0100] 最后，针对Loc_nv的Nnv个负值进行非负插值，具体步骤如下：
[0101] ①令t＝1:Nnv，并定义ψt为图像xa中以 为中心的大小为r×r的窗口，其中，表示图像xa在Loc_nv(t)处的负值像素。再对ψt中的非负像素按升序排序，记为ψ't；
[0102] ②如果 或者 即当 的95％置信区间的上边界 小于0或者 大于0但小于ψ't窗口中的最小像素值时，计算
得到ψ't窗口中符合要求的像素个数 其中， 表示正向取整操作；ρ的
作用是控制在非负插值中起作用的像素的百分比。再计算 表示
ψ't窗口中的像素集合，用 代替xfilter中对应位置上的值，得到平滑图像
[0103] ③如果不满足步骤②中的条件，则计算ψ't窗口中小于等于 的像素个数Nu，即满足 或者计算ψ't窗口中小于等于 且大于等于 的像素个数Nu，即满足 再计算 用 代
替xfilter中对应位置上的值，得到平滑图像
[0104] 步骤5：基于稀疏引导正则项约束的子目标函数，将步骤4得到的平滑图像 和步骤2估计的正则化参数τ作为该目标函数的暖启动，获得最终平滑且边缘信息丰富的高能X射线去模糊图像。
[0105] 步骤4得到的 和步骤2估计的正则化参数τ将会作为式(4)的暖启动，利用梯度投影稀疏重构算法(GPSR)(参照文献：Gradient Projection for Sparse Reconstruction:Application to Compressed Sensing and Other Inverse Problems)的改进方法——GPSR-BB算法对 进行非盲去模糊，得到最后恢复的去模糊图像。良好的起点对梯度投影法具有很大的作用。本发明使用步骤4得到的 和步骤2估计的正则化参数τ来初始化GPSR-BB，从而求得式(4)解附近的τ值。迭代的次数取决于τ值的近似值和解的近似值。利用这种暖启动技术，可以减少迭代的次数，并且有效地求解τ的一系列值。
[0106] GPSR方法的第一个关键步骤是将式(4)表示为一个二次规划方程，这可以通过将变量 分成正的和负的两部分来实现，形式上，则是引入向量并做替换：
[0107]
[0108] 对所有的i＝1,2,...,n，这些关系由 和 满足。其中，(·)+表示正向操作，即 因此有， 其中，1n＝[1,1,...,1]T。所
以，式(4)可以被重写为有界约束的二次规划方程(BCQP)：
[0109]
[0110] 虽然，这里的L2范数项不受u←u+s和v←v+s的影响。其中，s是移位向量。然而，这种变化增加了s这个变量。所以，式(20)可以用更标准的BCQP格式表示：
[0111]
[0112] 其中，
[0113]
[0114] GPSR-BB算法的具体步骤如下:
[0115] ①初始化：将步骤4得到的 作为初始值z(0)，选择标量参数αmin,αmax,α(0)∈[αmin,αmax]，并设置迭代次数n＝0。
[0116] ②计算：
[0117]
[0118] 其中， 的梯度导数。
[0119] ③线搜索：寻找标量λ(n)，使得F(z(n)+λ(n)δ(n))在区间λ(n)∈[0,1]最小化，并设置z(n+1)＝z(n)+λ(n)δ(n)。因为，F是二次方程式，所以线搜索参数λ(n)可以通过如下的解析表达式来计算：
[0120]
[0121] 当(δ(n))TBδ(n)＝0时，则令λ(n)＝1。引入λ(n)这个参数的作用是：消除了F在某些迭代中可能增大的可能性，它比GPSR的
基础方法(将λ(n)直接设置为1)的性能更好。
[0122] ④更新标量α：计算
[0123] γ(n)＝(δ(n))TBδ(n)； (25)
[0124] 如果γ(n)＝0，则令α(n+1)＝αmax，否则令
[0125] ⑤进行收敛性测试，并运用如下的终止准则进行判断，当近似解z(n+1)满足要求时停止迭代，否则令n←n+1并返回步骤②：
[0126]
[0127] 其中，εt为一个正的常数；tolA为设定的停止
阈值。
[0128] 本发明提出了一种高能X射线图像非盲去模糊方法。首先，在频域上运用Gibbs抽样动态构造马尔可夫链，进行快速非盲去模糊，能够提高MCMC算法的计算速度，降低计算的复杂度；然后，通过MCMC采样获得噪声标准差的统计量，可以作为迭代引导滤波步骤中迭代次数的判断依据；接着，基于MCMC统计样本95％的置信区间进行非负插值，剔除MCMC采样后图像中残留的负值像素，同时抑制图像的振铃效应；最后，结合经过引导滤波和非负插值处理后的平滑图像以及通过MCMC采样获得的正则化参数作为稀疏引导正则项约束的目标函数的暖启动进行非盲去模糊，获得最终的高能X射线去模糊图像。相较于其他方法，本发明可以加速高能X射线图像的去模糊进程，更好地去除图像中的噪声和振铃，在保持平滑度的同时获得更加精确的去模糊结果及其不确定度的统计量，而不会产生额外的噪声。
[0129] 本领域内的技术人员应明白，本
申请的实施例可提供为方法、系统、或
计算机程序产品。因此，本申请可采用完全
硬件实施例、完全
软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘
存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
[0130] 本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的
流程图和/或方
框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程
数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中
指定的功能的装置。
[0131] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0132] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0133] 以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和
变形，这些改进和变形也应视为本发明的保护范围。
