基于经验模态分解的激光散斑噪声抑制方法
1 引言
定量相位成像(QPI)[1-2]是一种新兴的无标记光学成像方法,可以实现反射样品三维形貌的可视化,测量透明和半透明样品的内部结构和折射率分布,其在透明生物样品上的应用是当前生物光学成像的研究热点之一。定量相位成像技术主要依赖全息[3]和干涉技术记录透明物体的相位信息,或是记录光场通过透明物体的衍射强度图,从记录的强度图中提取物体的相位。它可以与压缩感知[4]、深度学习[5-7]、自动对焦[8]、拓扑调制[9]或补偿技术[10]相结合,实现更高精度的测量结果。当用激光照射待测样品时,物体表面相对于激光波长量级而言相当粗糙,由于激光的高相干性,各个物点的散射光之间可能发生干涉,从而产生呈颗粒状无规则分布的散斑噪声。其存在会显著降低成像结果的分辨率、减少图像的信噪比,无法满足高精度相位重建的要求,因此有必要对激光散斑噪声进行抑制。
目前,常用的散斑噪声抑制的方法主要分成两类[11],一类是添加或替换硬件降低光源的相干性,如使用LED等低相干性的光源[12]等。Farrokhi[13]等人采用一动一静双散射片系统来获得高速散斑场照明,以抑制定量相位成像中的相干散斑。李煊[14]等人提出一种基于旋转双散射片的全息散斑噪声抑制方法,其光学信号的非相关周期大,为获得大量非相关全息图求均值去噪提供了依据。总体来说,这类硬件方法对光程差进行精确控制的要求更高,增加了搭建光路的复杂度和难度。另外一类则是通过图像处理手段实现散斑噪声抑制。Uzan等人[15]对位于图像块中心的有噪声像素采用非局部均值(NLM)滤波,对散斑噪声起到了较好的抑制效果。但该方法在去噪的同时平滑了部分边缘处的信息,会导致相位出现误差。刘吉等人[16]提出一种基于正余弦分解结合自适应全变分的去噪方法,该方法提高了图像的峰值信噪比,减少了去噪后相位图的波动性。但该方法中自适应参数的选取是图像去噪和保留边缘信息的关键,需要通过缩小自适应参数的取值范围来达到最优去噪的重建效果。Montresor等人[17]应用残差网络(DNCNN)对数字全息干涉测量的相位图进行去噪,具有非高斯统计和非平稳特性,并表现出空间相关长度和出色的去噪性能。但该方法需要使用大量的含散斑噪声的全息图数据集对多等级神经网络进行训练,实验条件复杂。牛瑞等人[18]提出利用二维高斯窗口对包裹相位进行散斑噪声抑制,不仅保留了包裹相位的跳变边缘,还提高了相位重建的精度。但该方法中二维高斯窗口大小以及阈值会影响相位重建精度,需要研究最佳窗口大小及阈值。吴育民等人[19]提出一种基于Canny算子改进各向异性扩散(P-M)方程的全息散斑抑制方法,在去噪的同时有利于保留图像的细节信息。但该方法的再现图像仍然存在轻微的二次干涉条纹,且图像相较于原始图像而言存在一定的模糊。这类图像处理手段虽然在定量相位成像中达到了一定的散斑抑制效果,但对于相位重建精度还有待提高。
针对以上问题,本文提出了一种结合经验模态分解(EMD)的Canny算子改进P-M方程的定量相位成像散斑噪声抑制方法。首先对所记录的全息图进行EMD得到一系列频率由高到低的固有模态函数(IMF)分量,仅用高频的IMF分量重构图像以增强图像细节信息。随后引入Canny算子对细节信息突出的重构图像进行边缘检测,从而更好地控制扩散去噪过程,提高了P-M方程抑制散斑噪声的性能,并通过实验验证了该散斑噪声抑制方法的可行性。
2 基本原理
2.1 结合EMD的Canny算子改进P-M方程法在全息散斑噪声抑制中的应用
在全息图中的散斑噪声是一种乘性噪声[20-21],其数学模型如
其中:
结合EMD的Canny算子改进P-M方程法首先将图像
(1)对图像
由于全息图的细节信息主要包含在高频分量
(2)对重构的图像
(3)根据边缘检测结果引导各向异性扩散方程的扩散程度,进行去噪,得到散斑抑制后的全息图
(4)随后对散斑抑制后的全息图
2.2 EMD
二维经验模态分解的目的是将一个二维信号
(1)找出
均值
(2)将
值得注意的是,
(3)令
(4)令
(5)重复上述过程
2.3 Canny算子改进P-M方程法
基于偏微分方程的图像去噪方法得到了广泛的应用,从线性均匀扩散最终发展到各向异性扩散。在各向异性扩散中,应用最为广泛的是Perona和Malik共同提出的Perona-Malik(P-M)方程[22],许多方法都是基于P-M方程发展起来的[23-24]。
然而传统的P-M方程法存在以下问题[25]:
(1)由于梯度算子
(2)梯度算子
为了克服上述问题,通过引入Canny边缘检测算子改进P-M方程从而克服梯度算子
其中:
在
使用Canny算子
(1)对图像
(2)按上、下、左、右4个方向求解图像
(3)令扩散系数
(4)由于Canny边缘检测结果
(5)散度算子
(6)图像
3 实验
3.1 实验装置及样本
本文搭建了一套基于Mach-Zehnder干涉方法的数字记录光路,如
图 1. 基于Mach-Zehnder干涉仪的数字记录光路原理图及实物图。(a)原理图;(b)实物图。
Fig. 1. Schematic diagram and physical diagram of digital recording optical path based on Mach-Zehnder interferometer.(a)Schematic diagram;(b)Physical diagram.
在如
图 2. 样品及全息图。(a)样品板;(b)样品A;(c)样品B;(d)样品C;(e)样品D;(f)样品E;(g)样品A对应的全息图A;(h)样品B对应的全息图B;(i)样品C对应的全息图C;(j)样品D对应的全息图D;(k)样品E对应的全息图E。
Fig. 2. Samples and holograms.(a)Sample plate;(b)Sample A;(c)Sample B;(d)Sample C;(e)Sample D;(f)Sample E;(g)Hologram A corresponding to sample A;(h)Hologram B corresponding to sample B;(i)Hologram C corresponding to sample C;(j)Hologram D corresponding to sample D;(k)Hologram E corresponding to sample E.
3.2 结合EMD的Canny算子改进P-M方程法的相位重建
首先,将全息图进行EMD,得到一系列从高频到低频的IMF分量,通过选取不同的IMF分量可以重构出不同效果的图像,
图 3. 选取不同IMF分量的重构结果。(a)选取 的重构图像;(b)选取 的重构图像;(c)选取 的重构图像;(d)选取 的重构图像。
Fig. 3. Reconstruction results with different IMF components selected.(a)Reconstructed image with ;(b)Reconstructed image with ;(c)Reconstructed image with ;(d)Reconstructed image with .
其次,用Canny算子对全息图进行边缘检测。
图 4. Canny算子边缘检测结果。(a)Canny算子对图3(a)的边缘检测结果;(b)Canny算子对图2(h)的边缘检测结果。
Fig. 4. Canny operator edge detection results.(a)Edge detection results of Canny operator for Fig.3(a);(b)Edge detection results of Canny operator for Fig.2(h).
我们用Canny算子改进P-M方程法对5张由最高频IMF分量重构的图像进行去噪,结果如
图 5. 经过所提出方法处理后的全息图。(a)全息图A;(b)全息图B;(c)全息图C;(d)全息图D;(e)全息图E。
Fig. 5. Holograms processed by the proposed method.(a)Hologram A;(b)Hologram B;(c)Hologram C;(d)Hologram D;(e)Hologram E.
最后,对处理后的全息图采用角谱衍射法[26]进行数值重建,从再现像中恢复出包裹相位信息。采用基于离散余弦的最小二乘法[27]对包裹相位进行展开,以此获得真实相位。针对所搭建的基于Mach-Zehnder干涉仪的数字记录光路,其相位误差的主要来源是物光与参考光之间引入的倾斜误差。采用最小二乘拟合法[28]进行相位补偿,最终得到的相位图如
图 6. 经过所提出方法处理后的重建相位图和未去噪全息图的重建相位图对比。(a)全息图A经过所提出方法处理后重建的相位图;(b)全息图B经过所提出方法处理后重建的相位图;(c)全息图C经过所提出方法处理后重建的相位图;(d)全息图D经过所提出方法处理后重建的相位图;(e)全息图E经过所提出方法处理后重建的相位图;(f)未去噪全息图A的重建相位图;(g)未去噪全息图B的重建相位图;(h)未去噪全息图C的重建相位图;(i)未去噪全息图D的重建相位图;(j)未去噪全息图E的重建相位图。
Fig. 6. Comparison of the reconstructed phase maps after processing by the proposed method and the reconstructed phase maps without denoising holograms.(a)Reconstructed phase map of hologram A after processing by the proposed method;(b)Reconstructed phase map of hologram B after processing by the proposed method;(c)Reconstructed phase map of the hologram C after processing by the proposed method;(d)Reconstructed phase map of hologram D after processing by the proposed method;(e)Reconstructed phase map of the hologram E after processing by the proposed method;(f)Reconstructed phase map of hologram A without denoising;(g)Reconstructed phase map of hologram B without denoising;(h)Reconstructed phase map of hologram C without denoising;(i)Reconstructed phase map of hologram D without denoising;(j)Reconstructed phase map of hologram E without denoising.
3.3 实验对比
为了验证所提出方法的有效性,将所提出的方法与基于Canny算子改进P-M方程法、均值滤波法以及中值滤波法进行比较。所记录的5种全息图经过基于Canny算子改进P-M方程法、均值滤波法和中值滤波法处理后,经过角谱衍射法数值重建[26]、基于离散余弦变换的最小二乘法相位展开[27]以及最小二乘拟合法[28]相位补偿之后得到的相位图如
图 7. 全息图A、B、C、D、E经过基于Canny算子改进P-M方程法、均值滤波法与中值滤波法处理后重建的相位图。(a)~(e)经过基于Canny算子改进P-M方程法处理后重建的相位图;(f)~(j)经过均值滤波法处理后重建的相位图;(k)~(o)经过中值滤波法处理后重建的相位图。
Fig. 7. Reconstructed phase maps of hologram A,B,C,D and E after processing by the improved P-M equation method based on the Canny operator,mean filter method and median filter method.(a)~(e)Reconstructed phase maps after processing by the improved P-M equation method based on Canny operator;(f)~(j)Reconstructed phase maps after processing by the mean filtering method;(k)~(o)Reconstructed phase maps after processing by the median filter method.
4 分析与讨论
4.1 降噪定量分析
为了对降噪后的重建相位图做定量分析,本文采用结构相似性(SSIM)[29]、边缘保持指数(EPI)[30]和散斑抑制指数(SSI)[31]作为定量分析的指标。
SSIM的计算公式如
其中:
EPI的计算公式如
其中:
SSI的计算公式如
其中:
5种全息图经过所提出方法、基于Canny算子改进P-M方程法、均值滤波法以及中值滤波法处理之后重建相位图对应的SSIM、EPI和SSI计算结果如
表 1. 不同散斑抑制方法的SSIM、EPI、SSI计算结果
Table 1. Calculation results of SSIM,EPI and SSI for different speckle suppression methods
|
从
以均值滤波对应的EPI为计算基准,中值滤波法相对于均值滤波法的EPI最大提高了15.027 8%,最小提高了1.487 6%,平均提高了6.630 7%;基于Canny算子改进P-M方程法相对于均值滤波法的EPI最大提高了23.961 8%,最小提高了1.999 7%,平均提高了12.449 1%;本文所提出方法相对于均值滤波法的EPI最大提高了24.606 4%,最小提高了3.392 5%,平均提高了14.386 1%。
以均值滤波的SSI为计算基准,中值滤波法相对于均值滤波法的SSI最大降低了8.942 5%,最小降低了0.072 3%,平均降低了2.918 3%;基于Canny算子改进P-M方程法相对于均值滤波法的SSI最大降低了9.548 8%,最小降低了2.546 9%,平均降低了5.089 3%;本文所提出方法相对于均值滤波法的SSI最大降低了12.125 5%,最小降低了5.763 5%,平均降低了8.129 9%。
4.2 相位截面曲线分析
为了验证降噪后定量相位重建的精确性,在如
图 8. 图6(a)、图7(a)、图7(f)和图7(k)所示相位图与原始相位图在y=711处的截面曲线对比。
Fig. 8. Phase maps shown in Fig.6(a),Fig.7(a),Fig.7(f)and Fig.7(k)compared with the cross-sectional curves of the original phase map at y=711.
图 9. 图6(b)、图7(b)、图7(g)和图7(l)所示相位图与原始相位图在x=415处的截面曲线对比。
Fig. 9. Phase maps shown in Fig.6(b),Fig.7(b),Fig.7(g)and Fig.7(l)are compared with the cross-sectional curves of the original phase map at x=415.
图 10. 图6(c)、图7(c)、图7(h)和图7(m)所示相位图与原始相位图在x=589处的截面曲线对比。
Fig. 10. Phase maps shown in Fig.6(c),Fig. 7(c),Fig.7(h)and Fig.7(m)compared with the cross-sectional curves of the original phase map at x=589.
图 11. 图6(d)、图7(d)、图7(i)和图7(n)所示相位图与原始相位图在y=959处的截面曲线对比。
Fig. 11. Phase maps shown in Fig.6(d),Fig.7(d),Fig.7(i)and Fig.7(n)compared with the cross-sectional curves of the original phase map at y=959.
图 12. 图6(e)、图7(e)、图7(j)和图7(o)所示相位图与原始相位图在x=1 355处的截面曲线对比。
Fig. 12. Phase maps shown in Fig.6(e),Fig.7(e),Fig.7(j)and Fig.7(o)compared with the cross-sectional curves of the original phase map at x=1 355.
从相位截面曲线变化来看,绿色虚线、蓝色点、橙色点线、紫色双点线与红色实线的变化趋势基本一致,但是绿色虚线与红色实线的偏差要比蓝色点、橙色点线、紫色双点线与红色实线的偏差要小,具体的偏差如
表 2. 绿色虚线、蓝色点、橙色点线和紫色双点线与红色实线的最大、最小、平均偏差
Table 2. Maximum,minimum and average deviation of the green dashed line,the blue point,the orange point line and the purple double point line from the red solid line
|
综合上述结果可以得出,在定量相位成像中,本文所提出的方法——结合EMD的Canny算子改进P-M方程法相对于基于Canny算子改进P-M方程法、均值滤波法和中值滤波法的SSIM、EPI与SSI计算结果以及截面曲线对比结果均更好。本文所提出方法具备更好的全息散斑噪声抑制效果,误差更小,相位重建精度更高。
5 结论
本文提出一种结合EMD的Canny算子改进P-M方程的散斑噪声抑制方法,该方法用来抑制在定量相位成像中高相干性光源产生的散斑噪声对相位重建精度的影响。通过EMD获得高频分量的IMF重构出细节突出的图像,使得Canny算子检测更加准确,从而更好地控制P-M方程的扩散去噪过程,在去噪的同时又保留高精度的相位信息。文中对本文所提出方法、基于Canny算子改进P-M方程法、均值滤波法和中值滤波法降噪后的重建相位图用SSIM、EPI和SSI进行了定量分析。其中,以均值滤波为参考,本文所提方法去噪后的重建相位图像的SSIM提高了12.900 0%,EPI提高了14.386 1%,SSI降低了8.129 9%。通过对不同方法降噪处理后重建的相位图像进行截面曲线分析,经对比发现,本文所提方法与原始相位偏差更小。实验结果表明,本文所提方法在性能上优于其他3种降噪方法,有较好的散斑抑制效果,相位重建精度更高。
[1] POPESCUG. Quantitative Phase Imaging of Cells and Tissues [M]. New York: McGraw-Hill, 2011: 121-138.
[14] 李煊. 基于旋转双散射片的数字全息显微成像散斑抑制研究 [D]. 西安: 西安理工大学, 2020. 10.30919/esee8c722
[16] 刘吉, 黄晓慧, 武锦辉, 等. 基于正余弦分解的自适应全变分散斑去噪方法[J]. 中国激光, 2020, 47(10): 1004004.
[18] 牛瑞, 田爱玲, 王大森, 等. 数字全息测量系统的散斑噪声抑制[J]. 激光与光电子学进展, 2022, 59(16): 1609002.
[19] 吴育民, 段海燕, 文永富, 等. 再现图像细节抑制散斑噪声技术研究[J]. 影像科学与光化学, 2018, 36(2): 187-191.
[21] 王灿进, 石宁宁, 孙涛. 同态非局部滤波在激光主动成像散斑抑制中的应用研究[J]. 液晶与显示, 2016, 31(2): 193-200.
[25] 陈一虎. P-M扩散方程图像去噪方法分析[J]. 宝鸡文理学院学报(自然科学版), 2010, 30(4): 14-17.
CHEN Y H. An analysis of the method of P-M diffusion equation for image denoising[J]. Journal of Baoji University of Arts and Sciences (Natural Science Edition), 2010, 30(4): 14-17.
[27] 彭景, 张蓉竹. 干涉检测中常用解包裹算法噪声处理能力比较[J]. 光学与光电技术, 2019, 17(4): 77-83.
[28] 马树军, 刘炜华, 周鹏飞. 一种离轴数字全息显微相位自动补偿方法[J]. 东北大学学报(自然科学版), 2019, 40(6): 847-851.
MA S J, LIU W H, ZHOU P F. Off-Axis digital holographic microscopic phase automatic compensation method[J]. Journal of Northeastern University (Natural Science), 2019, 40(6): 847-851.
[30] 金鑫, 张景雄. 自适应加权中值滤波的InSAR干涉图去噪方法[J]. 测绘地理信息, 2016, 41(3): 12-15.
JIN X, ZHANG J X. A self-adaptive weighted-median filter for InSAR interferograms[J]. Journal of Geomatics, 2016, 41(3): 12-15.
[31] 汤春明, 张洪科, 于翔, 等. 反射式离轴数字全息再现像中的噪声抑制[J]. 半导体光电, 2014, 35(4): 745-748.
Article Outline
詹晓江, 甘楚立, 丁毅, 胡轶, 许彬, 习江涛, 邓定南. 基于经验模态分解的激光散斑噪声抑制方法[J]. 液晶与显示, 2023, 38(4): 495. Xiao-jiang ZHAN, Chu-li GAN, Yi DING, Yi HU, Bin XU, Jiang-tao XI, Ding-nan DENG. Laser speckle noise suppression method based on empirical mode decomposition[J]. Chinese Journal of Liquid Crystals and Displays, 2023, 38(4): 495.