DOI： 10.11817/j.issn.1672-7207.2021.02.013

基于非局部Bayesian的单掩膜编码衍射成像算法

司菁菁^{1, 2}，党文伟¹，史姗姗¹，程银波³

(1. 燕山大学 信息科学与工程学院，河北 秦皇岛，066004；

2. 河北省信息传输与信号处理重点实验室，河北 秦皇岛，066004；

3. 河北农业大学 海洋学院，河北 秦皇岛，066003)

摘要：为了降低编码衍射成像系统的硬件复杂度，研究基于单个编码衍射图案(coded diffraction pattern， CDP)的成像算法。以图像本身具有的非局部区域自相似特性作为相位恢复中的信号约束，克服单个CDP携带信息过少的不足，提出一种基于非局部Bayesian估计与图像自相似性的相位恢复与图像重建算法。研究结果表明：当仅利用单个CDP重建图像时，本文算法获得的平均峰值信噪比分别比WF(Wirtinger flow)，DOLPHIn(dictionary learning for phase retrieval)和ConPR(constrained phase retrieval)算法提高了15.00，5.00和0.23 dB；且其重建图像能够更好地保留细节信息，具有较高的主观视觉质量。此外，本文算法的运行时间约为DOLPHIn算法的1/4、ConPR算法的一半。

关键词：相位恢复；编码衍射图案；图像自相似性；Bayesian估计

中图分类号：TN919.8                    文献标志码：A

文章编号：1672-7207（2021）02-0450-08

Coded diffraction imaging algorithm with a single mask based on non-local Bayesian

SI　Jingjing^{1, 2}, DANG　Wenwei¹, SHI　Shanshan¹, CHENG　Yinbo³

(1. School of Information Science and Engineering, Yanshan University, Qinhuangdao 066004, China;

2. Hebei Key Laboratory of Information Transmission and Signal Processing, Qinhuangdao 066004, China;

3. Ocean College, Hebei Agricultural University, Qinhuangdao 066003, China)

Abstract: To reduce the hardware complexity of coded diffraction imaging system, the coded diffraction imaging algorithm was studied based on a single coded diffraction pattern(CDP). Taking the non-local self-similarity characteristics of the image as the regularization in the process of phase retrieval, the lack of information carried by a single CDP was offset. A new phase recovery and image recovery algorithm was proposed based on the non-local Bayesian estimation and the self-similarity of the image. The results show that when only one CDP is involved, the peak signal to noise ratio achieved by our algorithm is in average 15.00, 5.00 and 0.23 dB higher than those achieved by WF(Wirtinger flow), DOLPHIn(dictionary learning for phase retrieval) and ConPR(constrained phase retrieval), respectively. Images reconstructed by our algorithm maintain more detailed information, and present higher subjective visual quality. Moreover, the runtime of our algorithm is approximately a quarter of that of DOLPHIn and half of that of ConPR.

Key words: phase retrieval; coded diffraction pattern; image self-similarity; Bayesian estimation

利用物体在某种射线照射下的衍射光传播到探测器上的数据对物体进行高分辨率成像，是材料学和生物医学等诸多领域的一种重要研究手段。因为目前探测器能直接记录的通常仅是衍射场的强度信息，所以，在大多数情况下，衍射场的相位信息会在数据记录过程中丢失^[1]。然而，这些丢失的相位通常包含比强度更多的物体结构信息，对于物体成像与检测至关重要。因此，研究利用探测器测量到的衍射场强度信息恢复相位信息从而重建信号的相位恢复(phase retrieval)技术具有十分重要的意义。

2015年，CANDES等^[2]设计了编码衍射成像系统(coded diffraction imaging system)。与传统衍射成像系统不同，编码衍射成像系统在待成像物体之后加入了随机相位板(掩膜)，探测器记录的是编码衍射图案(coded diffraction patterns, CDP)。与传统的强度测量值相比，CDP中包含了更丰富的图像信息。此后，基于CDP的相位恢复与图像重构逐渐成为相位恢复领域的研究热点。CANDES等^[3]提出了采用谱估计初始化的WF(Wirtinger flow)算法，利用梯度下降迭代出最优解。CHEN等^[4]提出了TWF(truncated Wirtinger flow)算法，引入截断机制，促进梯度下降向更准确的方向进行。WANG等^[5]提出了结合截断谱估计初始化和截断振幅的TAF(truncated amplitude flow)算法。这些算法未利用图像内在的先验信息，需要较多数量的CDP才能实现图像的清晰重构，且算法的抗噪性能较差。实际上，CDP数量越少(即需要加入的掩膜数量越少)，成像系统的硬件越简单，越符合实现应用的需要。然而，CDP数量越少，包含的图像信息越少，重建原始图像的难度越大。GROSS等^[6]在编码衍射成像系统中引入PhaseLift技术，明显降低了重建图像所需的CDP的数量，但算法的运行时间显著增加。TILLMANN等^[7]利用自适应合成字典下的稀疏表示模型构建相位恢复优化问题中的正则项，提出了DOLPHIn(dictionary learning for phase retrieval)算法。CHANG等^[8-9]在优化模型中利用全变差(total variation, TV)正则化，进而设计了基于正交字典学习的相位恢复算法。SHI等^[10-11]提出了利用梯度稀疏与结构稀疏的双稀疏正则化模型(double sparse regularization, DSR)和一种结合交替投影与正则化的约束相位恢复(constrained phase retrieval, ConPR)算法。这些研究表明，在编码衍射成像算法中恰当地引入图像特有的先验信息，能够利用较少数量的CDP实现图像的较高质量重建。

自然图像本身具有非局部区域自相似性特征，其非局部区域纹理和结构会多次重复出现。若能在解决图像反问题时充分利用这一特性，则能更好地保持图像的纹理和细节信息^[12-14]。LEBRUN等^[15]基于图像非局部自相关性，提出了非局部Bayesian(non-local Bayesian, NLB)算法，根据Bayesian最大后验概率估计，建立目标图像块与其相似图像块集合的约束关系。本文基于图像自相似性与非局部Bayesian估计^[15-16]的约束关系引入编码衍射成像系统中，研究基于单个CDP的编码衍射成像算法，利用自然图像本身具有的非局部区域自相似特性构造正则项，进而构建整体的相位恢复优化问题，克服单个CDP携带信息过少的不足，提出一种基于非局部Bayesian估计与图像自相似性的相位恢复算法(phase recovery algorithm based on non-local Bayesian estimation and image self-similarity, prNLBEIS)。

1  相位恢复与编码衍射成像

相位恢复是非线性压缩感知理论的一个分支，它的实质是利用探测器测量到的衍射场强度信息恢复相位信息，从而重建信号。在没有噪声干扰的理想状态下，相位恢复问题的幅值测量可以表示为

(1)

若在测量过程中受到了加性噪声的干扰，则幅度测量可以表示为

(2)

式中：为待重建的原始信号；为幅值测量向量；为线性测量矩阵；为加性噪声。

在编码衍射成像系统中，测量矩阵为

(3)

式中：为傅里叶变换矩阵；M_L为掩膜矩阵；为掩膜矩阵的数量，即CDP的数量。越大，CDP中携带的待重建图像的信息越多，但相应的编码衍射成像系统的实现成本越高。本文研究基于单个CDP的相位恢复算法，即令，此时。

图像本身含有大量的先验信息，为了能够在相位恢复中充分利用图像自身的特性，可以将图像先验特性作为优化问题的约束，利用先验信息构造优化问题中的正则项。正则化的相位恢复模型可以表示为

(4)

式中：为数据保真项；为包含图像先验信息的正则项；为正则项的控制参数。

2  基于非局部Bayesian与图像自相似性的相位恢复算法

2.1　基于非局部Bayesian与图像自相似性约束的图像恢复

将图像的非局部自相似特性作为图像的先验信息，根据Bayesian最大后验概率估计，建立目标图像块与其相似图像块集合的约束关系，并在此基础上实现图像恢复。

令x表示原始图像，表示受高斯白噪声污染后的退化图像。图像恢复问题可以表示为

(5)

由Bayesian定理可知：

(6)

从而将图像恢复问题转换成

(7)

其中：

(8)

式中：为高斯白噪声的标准差；为归一化常数。

为了求解如式(7)所示的图像恢复问题，需要利用原始图像x的先验信息。本文利用自然图像的非局部区域自相似特性，从而在重建图像中尽可能地保持原始图像中的纹理和细节信息。

将原始图像x划分成大小相等的若干图像块，并用r表示其中一块。如式(7)所示的图像x的恢复问题可转换为若干个如式(9)所示的恢复图像块r的子问题：

(9)

其中，

(10)

设与r相似的图像块服从高斯分布，即

(11)

式中：为与r相似的图像块集合的均值；为与r相似的图像块集合的协方差矩阵。在此基础上，式(9)可以转换为

(12)

式(12)中的和需要根据原始图像计算，在实际图像恢复过程中并不能直接得到。

为了解决此问题，对退化图像进行等尺寸分块，并利用块匹配方法为每个图像块搜索出相似块集合。根据加性高斯白噪声的退化模型，利用退化块的相似块集合均值和协方差矩阵分别近似计算出原始图像块r对应的与：

(13)

(14)

这样，可将如式(12)所示的优化问题转换为

(15)

将式(15)中的优化函数对r求导，并令其等于0，得到

(16)

整理化简可得

(17)

实现了由退化图像块对原始图像块的估计。

利用式(17)对每个退化图像块进行处理之后，将会得到每个图像块r的初步重建结果。若利用这些初步重建块重新进行块匹配，则能更准确地搜索出与相似的图像块，且在其基础上计算出的相似块集合均值与协方差矩阵与原始图像块r对应的与更接近。据此，本文采用二次估计的方法，基于非局部Bayesian估计与图像自相似性约束，实现图像块的重构。具体方法可总结为：

1) 将退化图像划分成大小相等的图像块并进行块匹配，为每个退化图像块的相似块集合计算均值和协方差矩阵，进而利用式(17)计算每个图像块的初步重建块。

2) 根据初步重建块重新进行块匹配，并为每个初步重建块的相似块集合计算均值和协方差矩阵，将和代入式(17)，得到式(18)，以计算每个图像块的二次重建结果。

(18)

3) 由所有的二次重建图像块组合重建图像。

2.2　基于非局部Bayesian与图像自相似性约束的单CDP成像算法

将基于非局部Bayesian估计与图像自相似性的约束关系引入编码衍射成像系统，实现单编码衍射图案下的相位恢复。

令表示对图像x进行的分块操作。设图像块的个数为n，将第i个图像块表示为，。根据式(15)，将基于非局部Bayesian估计与图像自相似性约束的优化函数表示为

(19)

式中：z为近似图像；和分别为z中第i个图像块的相似集均值与协方差矩阵。令，，则第i个图像块对应的优化函数可以表示为

(20)

将式(19)作为正则项代入式(4)，即可得到基于非局部Bayesian与图像自相似性约束的相位恢复模型：

(21)

式中：。

本文利用向前向后分裂(Forward-Backward Splitting, FBS)算法^[17]求解该优化问题。对于，向前梯度下降(forward gradient decent)

(22)

其中：

(23)

为矩阵的共轭转置矩阵；表示对应元素相乘。对于，向后梯度下降(backward gradient descent)，计算近端算子(proximal operator)为

(24)

根据式(19)，式(24)可以表示为

(25)

从而可将分解为n个子近端算子，其中，第i个可以表示为

(26)

对于如式(26)所示的优化问题，需要联合优化和。

基于以上分析，采用交替最小化的方法重建图像，具体如下：

1) 固定，优化：根据图像块组合出图像x，根据式(22)更新z，分块后得。

2) 固定，优化：求解形如式(12)所示的优化问题。根据进行块匹配，并为各图像块的相似块集合计算均值与协方差矩阵，从而根据式(27)获得非局部Bayesian的初步估计解：

(27)

进一步，根据重新进行块匹配，并计算各图像块相似集合的均值与协方差矩阵，利用式(28)获得第二次非局部Bayesian估计的结果：

  (28)

采用上述方法交替更新和，直到满足停止条件为止。将最终得到的经过块组合后，得到重建图像，从而实现如式(21)所示优化问题的求解，完成相位恢复与图像重建过程。

本文设计的prNLBEIS算法的具体实现过程如下。

输入：幅度测量向量，测量矩阵，迭代次数，误差容限，步长，噪声均方差修正系数。

输出：重构图像。

步骤如下。

初始化：

For (；t≤T；t++) do

{ 1)

2) 对进行一级二维小波变换，得小波系数矩阵，利用对角线子带HH中系数的中值，估计：

3) 对分块，得到各图像块，。

4) 根据进行块匹配，计算各图像块相似集合的均值与协方差矩阵，并在此基础上计算图像块的初步估计结果：

；

5) 根据图像块的初步估计结果进行块匹配，计算各图像块相似集合的均值与协方差矩阵，实现图像块的二次估计：

；

6) 由图像块组合出图像；

7) 若，则跳出循环。}

3  实验结果与分析

选用512像素×512像素的“Barbara”和“Lena”以及256像素×256像素的“Cameraman”和“peppers”4幅标准灰度图像作为实验对象，利用Matlab编写实验程序，验证本文设计的prNLBEIS算法的性能，并与现有的WF算法、DOLPHIn算法以及当前重构性能较高的ConPR算法进行比较。WF，DOLPHIn和ConPR这3种算法的Matlab代码均由原作者个人网站下载。进行基于单CDP的相位恢复与图像重构的实验，即令L=1。考虑高斯测量噪声的影响，根据测量值向量，重建图像。采用峰值信噪比(peak signal to noise ratio, PSNR)和结构相似性(structural similarity index measurement, SSIM)作为评价算法重建性能的客观指标。PSNR的计算公式为

(29)

式中：为峰值信噪比；为原始图像；为重建图像；为图像的像素；为像素的最大值。图像重构质量越高，越大。

SSIM的计算公式为

(30)

式中：为结构相似性指标；和分别为原始图像和重建图像的均值；和分别为和的标准差；为与的协方差；常数，，用以避免分母为0。图像重构质量越高，越大。

实验中，选用八元随机掩膜，即令掩膜矩阵M中有4/5的元素为{}，其余1/5的元素为{}。在本文prNLBEIS算法中，以高斯随机矩阵初始化。根据经验，实验取步长u=0.28，迭代次数=30，误差容限=0.000 1，噪声均方误差修正系数=0.004 9。

表1　不同信噪比下重建图像峰值信噪比比较

Table 1　Comparison of PSNR obtained under different SNR    峰值信噪比/dB

表2　不同信噪比下重建图像结构相似性指标比较

Table 2　Comparison of SSIM obtained under different SNR

在不同信噪比(signal to noise ratio, SNR)下，比较prNLBEIS，WF，DOLPHIn和ConPR这4种算法重建图像的客观指标。表1和表2所示为这4种算法在SNR分别为5，10，15和20 dB时重建图像的峰值信噪比和结构相似性指标。由表1和表2可见，在不同噪声水平下，prNLBEIS算法重建图像的峰值信噪比和结构相似性指标均明显比WF算法和DOLPHIn算法的高；在大多数情况下，prNLBEIS算法的峰值信噪比和结构相似性指标也比ConPR算法的高。与WF，DOLPHIn和ConPR算法相比，prNLBEIS算法获得的峰值信噪比平均值分别提高了15.00，5.00和0.23 dB。

由于WF算法在单CDP条件下不能有效地重构图像，因此，从主观视觉效果上比较DOLPHIn，ConPR算法与本文prNLBEIS算法的重建图像。当信噪比为10 dB时，DOLPHIn，ConPR算法与本文prNLBEIS算法的重构图像如图1所示。由图1可以看出，DOLPHIn算法的重构图像中存在着明显的块效应；ConPR算法的重建图像中缺失了大量的细节信息；prNLBEIS算法在重建图像中保持了更多的细节信息。

图1　SNR为10时重构的Barbara和Cameraman图像

Fig. 1　Reconstructed Barbara and Cameraman under SNR of 10

在信噪比为10 dB和20 dB下进行图像重建实验，比较4种算法的运行速度。以信噪比为10 dB时DOLPHIn算法重建256像素×256像素Cameraman图像的运行时间作为标准，计算其他算法的运行时间与该时间的比值，结果如表3所示。由表3可见，WF算法的运行时间最短；prNLBEIS算法的运行时间约为DOLPHIn算法运行时间的1/4；prNLBEIS算法的运行时间约为ConPR算法运行时间的一半。

表3　SNR为10 dB和20 dB时算法相对运行时间的比较

Table 3　Comparison of relative running time when SNR is 10 dB or 20 dB

综合以上实验结果可见，虽然WF算法的运行时间最短，但是其在单CDP的条件下不能有效地重构图像。prNLBEIS算法重建图像的客观指标与主观质量均优于DOLPHIn与ConPR算法，且prNLBEIS算法的运行时间明显低于DOLPHIn与ConPR算法的运行时间。

4  结论

1) 在编码衍射成像系统中，重建算法需要的CDP数量越少，成像系统的硬件越简单，越符合实现应用的需要。

2) 基于单个CDP的编码衍射成像，提出一种基于非局部Bayesian估计与图像自相似性的相位恢复prNLBEIS算法。

3) 与现有算法相比，prNLBEIS算法仅需单个CDP即能实现图像的较高质量重构，且具有较快的运行速度。

参考文献：

[1] SHECHTMAN Y, ELDAR Y C, COHEN O, et al. Phase retrieval with application to optical imaging: a contemporary overview[J]. IEEE Signal Processing Magazine, 2015, 32(3): 87-109.

[2] CANDES E J, LI Xiaodong, SOLTANOLKOTABI M. Phase retrieval from coded diffraction patterns[J]. Applied and Computational Harmonic Analysis, 2015, 39(2): 277-299.

[3] CANDES E J, LI X, SOLTANOLKOTABI M. Phase retrieval via Wirtinger flow: theory and algorithms [J]. IEEE Transaction on Information Theory, 2015, 61(4): 1985-2007.

[4] CHEN Yuxin, CANDES E J. Solving random quadratic systems of equations is nearly as easy as solving linear systems[J]. Communications on Pure and Applied Mathematics, 2017, 70(5): 822-883.

[5] WANG Gang, GIANNAKIS G B, ELDAR Y C. Solving systems of random quadratic equations via truncated amplitude flow[J]. IEEE Transaction on Information Theory, 2016, 64(2): 773-794.

[6] GROSS D, KRAHMER F, KUENG R. Improved recovery guarantees for phase retrieval from coded diffraction patterns[J]. Applied and Computational Harmonic Analysis, 2017, 42(1): 37-64.

[7] TILLMANN A M, ELDAR Y C, MAIRAL J. DOLPHIn: dictionary learning for phase retrieval[J]. IEEE Transactions on Signal Processing, 2016, 64(24): 6485-6500.

[8] CHANG Huibin, LOU Yifei, NG M K, et al. Phase retrieval from incomplete magnitude information via total variation regularization[J]. SIAM Journal on Scientific Computing, 2016, 38(6): A3672-A3695.

[9] CHANG Huibin, MARCHESINI S. Denoising Poisson phaseless measurements via orthogonal dictionary learning[J]. Optics Express, 2018, 26(16): 19773-19796.

[10] SHI Baoshun, LIAN Qiusheng, CHEN Shuzhen, et al. Coded diffraction imaging via double sparse regularization model[J]. Digital Signal Processing, 2018, 79: 23-33.

[11] SHI Baoshun, LIAN Qiusheng, HUANG Xin, et al. Constrained phase retrieval: when alternating projection meets regularization[J]. Journal of the Optical Society of America B: Optical Physics, 2018, 35(6): 1271-1281.

[12] CHANG K, DING P L K, LI Baoxin. Single image super-resolution using collaborative representation and non-local self-similarity[J]. Signal Processing, 2018, 149: 49-61.

[13] CHEN Honggang, HE Xiaohai, QING Linbo, et al. Single image super-resolution via adaptive transform-based nonlocal self-similarity modeling and learning-based gradient regularization[J]. IEEE Transactions on Multimedia, 2017, 19(8): 1702-1717.

[14] DE BORTOLI V, DESOLNEUX A, GALERNE B, et al. Patch redundancy in images: a statistical testing framework and some applications[J]. SIAM Journal on Imaging Sciences, 2019, 12(2): 893-926.

[15] LEBRUN M, BUADES A, MOREL J M. A nonlocal Bayesian image denoising algorithm[J]. SIAM Journal on Imaging Sciences, 2013, 6(3): 1665-1688.

[16] SHI Guiling, LIM C Y, MAITI T. Bayesian model selection for generalized linear models using non-local priors[J]. Computational Statistics & Data Analysis, 2019, 133: 285-296.

[17] BELLO CRUZ J Y, NGHIA T T A. On the convergence of the forward-backward splitting method with linesearches[J]. Optimization Methods and Software, 2016, 31(6): 1209-1238.

(编辑  赵俊)

收稿日期： 2020 -07 -09; 修回日期： 2020 -10 -07

基金项目(Foundation item)：国家自然科学基金资助项目(61701429)；河北省自然科学基金资助项目(F2018203137) (Project(61701429) supported by the National Natural Science Foundation of China; Project(F2018203137) supported by the Natural Science Foundation of Hebei Province)

通信作者：司菁菁，博士，副教授，从事多媒体信号处理研究；E-mail: sjj@ysu.edu.cn

引用格式： 司菁菁, 党文伟, 史姗姗, 等. 基于非局部Bayesian的单掩膜编码衍射成像算法[J]. 中南大学学报(自然科学版), 2021, 52(2): 450-457.

Citation: SI Jingjing, DANG Wenwei, SHI Shanshan, et al. Coded diffraction imaging algorithm with a single mask based on non-local Bayesian[J]. Journal of Central South University(Science and Technology), 2021, 52(2): 450-457.