光电工程  2019, Vol. 46 Issue (6): 1      DOI: 10.12086/oee.2019.180531     
结合多尺度分解和梯度绝对值算子的显微图像清晰度评价方法
崔光茫1,2,3 , 张克奇1 , 毛磊1 , 徐之海2 , 冯华君2     
1. 宁波永新光学股份有限公司,浙江 宁波 315040;
2. 浙江大学现代光学仪器国家重点实验室,浙江 杭州 310027;
3. 杭州电子科技大学电子信息学院,浙江 杭州 310018
摘要:针对显微图像自动对焦和成像系统质量评价问题,结合多尺度分解工具和梯度绝对值算子设计,提出了一种显微图像清晰度评价算法。采用非下采样剪切波分解,对输入的显微图像进行多尺度、多方向变换,得到一幅低频子带图像和若干幅高频子带图像。结合抗噪阈值设置,计算各子带图像的梯度绝对值算子和,利用图像清晰度变化对于低频和高频子带系数影响的差异,将高低频梯度绝对值算子的比值作为最终的显微图像清晰度评价数值。开展了仿真论证实验和实拍论证实验,实验结果表明,所提出的清晰度评价算法具有较好的单调性和抗噪特性,和几种经典的评价算法相比,本文方法得到的评价结果在灵敏度、稳定性和鲁棒性方面表现更为优异,具有很好的实际应用价值。
关键词显微图像    清晰度评价    非下采样剪切波变换    梯度绝对值算子    
Micro-image definition evaluation using multi-scale decomposition and gradient absolute value
Cui Guangmang1,2,3, Zhang Keqi1, Mao Lei1, Xu Zhihai2, Feng Huajun2     
1. Ningbo Yongxin Optics Co., Ltd., Ningbo, Zhejiang 315040, China;
2. State Key Lab of Modern Optical Instrumentation, Zhejiang University, Hangzhou, Zhejiang 310027, China;
3. School of Electronics and Information, Hangzhou Dianzi University, Hangzhou, Zhejiang 310018, China
Abstract: Aimed at the problem of automatic focus and image system quality evaluation in microscopy imaging, a micro-image definition evaluation method is presented by combining multi-scale decomposition tools and absolute gradient operators. The multiscale and multidirectional non-subsampled Shearlet transform is utilized to decompose the input micro image into a low frequency sub-band image and a number of high frequency sub-band images. Combined with the anti-noise threshold setting, the gradient absolute sum values of each sub-band image were calculated. By using the different effects of image sharpness on the low-frequency and high-frequency sub-band coefficients, the ratio of the high-frequency to low-frequency gradient absolute value operator was taken as the final evaluation value of the microscopic image sharpness. The simulation experiment and actual experiments were carried out and the experimental results illustrated that the proposed approach has good monotonicity and anti-noise characteristics. Compared with other classic evaluation algorithms, the presented method obtained superior performance on sensitivity, stability and robustness. It has very good practical application values.
Keywords: micro-images    image definition evaluation    non-subsampled Shearlet transform (NSST)    gradient absolute value    

1 引言

光学显微镜作为观察微观世界的重要仪器设备,已被广泛地应用在医疗健康、生物检测、工业生产等相关领域。由于显微镜成像系统景深非常小,只有在微小深度变化范围内才能采集得到清晰图像。显微图像清晰度的评价与判定,直接影响显微镜自动对焦的精度,也成为衡量显微系统成像质量的重要指标[1-3]。随着多媒体技术和数字图像的发展,显微仪器设备对于自动化的要求也逐渐提高,基于图像处理的显微图像清晰度评价算法得到越来越多的关注,这对于实现快速准确的显微自动对焦和成像系统性能评价具有重要意义[4-6]

按照通用的准则,图像清晰度评价方法一般可以分为主观评价和客观评价两大类。主观评价方法是最可靠、准确和全面的图像清晰度评价方法,人是成像链路的最终接受者,人眼视觉系统的主观感受才是评价图像质量好坏的标准。但是主观评价耗费大量的人力和时间,往往稳定性不够,不利于实时在线的评价和系统自动化的反馈处理。因此,客观评价方法成为实际应用中人们研究的重点。

近年来,研究者们在显微图像客观清晰度评价方法领域开展了较为深入的研究,从算法设计思路来看,可以分为基于梯度信息的方法、基于频域分解的方法以及基于熵函数的方法。对于基于梯度信息的方法,其大多数针对目标的细节边缘进行检测计算,利用边缘梯度函数统计值来判定图像的清晰程度,是应用较为广泛的算法。目前常用的梯度函数有Tenengrad函数[7]、方差(Variance)函数[8]、灰度平均梯度法(gray mean gradient, GMG)[9]、能量梯度(energy of gradient, EOG)函数[10]、Brenner函数、拉普拉斯算子和(Laplacian sum, LS)[11]、Robert函数等。而频谱函数评价方法是将图像变换到频域,认为图像的边缘与细节信息为频域高频分量,通过利用频域信息计算图像的清晰度评价值,常用的频域函数有傅里叶变换法、离散余弦变换(discrete cosine transform, DCT)函数和小波(Wavelet)变换。而利用熵函数的评价方法,通过统计图像灰度信息丰富程度来反映出图像清晰程度,算法计算简单,但往往容易受到噪声的影响,灵敏度不高[11]。近年来,多尺度分解工具被应用于图像清晰度评价,研究者们提出了一系列算法实施框架,取得了较好的结果,但算法在抗噪性和灵敏度等方面的性能表现也都有各自的局限性[5, 12-13]

显微图像获取过程中,容易受到成像链路噪声的影响,传统的清晰度评价方法往往会出现多个伪峰值,在算法灵敏度、无偏性、单峰性等方面有各自的局限性,给数字显微图像清晰度评价提出了新的挑战。本文结合多尺度分解工具和梯度绝对值算子设计,提出了一种显微图像清晰度评价算法。利用非下采样剪切波变换(non-subsampled Shearlet transform, NSST)将图像边缘信息分解到不同尺度的细节图层,通过抗噪阈值的设置,计算高频和低频的梯度绝对值能量算子,并利用高低频能量信息的变化特点,实现了有效的显微清晰度评价。设计了仿真实验和实拍实验,通过和其他几种评价算法进行对比论证,实验结果表明,本文算法具有更优的评价性能表现。

2 结合多尺度分解和梯度绝对值算子的清晰度评价

显微图像由于采集过程中视场范围较小,图像内容信息多少不确定,图像背景往往会对清晰度评判造成扰动,同时图像获取过程存在成像链路噪声,也会给清晰度评价方法带来新的挑战。频域多尺度分解工具是一种重要的图像分析处理方法,通过对图像进行分解,统计不同尺度分解图层的信息特征变化,进而建立与主观评价相关性较好的图像清晰度评价算法。本文引入非下采样剪切波变换(NSST),得到显微图像不同尺度的子带系数,并构建了梯度绝对值能量算子,利用低频和高频信息对图像清晰度变化的敏感程度差异,实现了高精度的显微图像清晰度评价算法。

2.1 非下采样Shearlet变换(NSST)

非下采样剪切波变换(NSST)是近年来提出的非常有效的图像多尺度分解工具,其数学结构简单,相比于曲波(Curvelet)和轮廓波(Contourlet)变换,具有抛物线尺度、方向敏感性更强及最优稀疏等特点[14]。同时,NSST避免了下采样过程对于图像信息的丢失,具有平移不变性和稳定性的优点,能够更好地表达出图像的轮廓、边缘、纹理等细节信息,适合图像特征的提取,可以为清晰度评价算法提供更多的判断信息。

根据合成小波理论,当分解维度n=2时,剪切波仿射系统函数可以表示为

${T_{{\boldsymbol{AB}}}}(y) = \{ {y_{j, l, k}} = {\left| {\det {\boldsymbol{A}}} \right|^{j/2}}y({{\boldsymbol{B}}^l}{{\boldsymbol{A}}^j}x - k)\}, $ (1)

其中:$j, l \in Z$$k \in {Z^2}$$y \in {L^2}({R^2})$L表示可积空间,AB是二阶可逆矩阵,${{\boldsymbol{A}}^j}$表示尺度变换矩阵,${{\boldsymbol{B}}^l}$表示几何变换矩阵,$\left| {\det {\boldsymbol{B}}} \right| = 1$。当${T_{{\boldsymbol{AB}}}}(y)$满足紧支撑(Parseval)框架时,函数系统里的元素称为合成小波,表示为

${\sum\limits_{j, l, k} {\left| {\left\langle {f, {y_{j, l, k}}} \right\rangle } \right|} ^2} = {\left\| f \right\|^2}, \forall f \in {L^2}({R^2})。$ (2)

${{\bf{A}}_a} = \left[ \begin{gathered} a 0 \\ 0 \sqrt a \\ \end{gathered} \right]$${{\bf{B}}_s} = \left[ \begin{gathered} 1 s \\ s 1 \\ \end{gathered} \right]$时,多分辨率多方向

性的Shearlet变化系统表示为

$\left\{ {{y_{ask}}(x) = {a^{ - \frac{3}{4}}}y({\bf{A}}_a^{ - 1}{\bf{B}}_s^{ - 1}x - k), a \in {R^ + }, s \in R, k \in {R^2}} \right\}, $ (3)

其中Shearlet变换函数中有3个参数变量:$a, s, k$$a$为图像尺度大小参数,s为剪切方向参数,k为平移量参数。通过3个参数的设置,可以表征不同方向、不同位置、不同尺度上的图像细节信息。

NSST是Shearlet变换的无下采样变换方式,算法多尺度分解流程图,如图 1所示,其主要处理步骤包括:

图 1 非下采样剪切波分解流程图 Fig. 1 Flowchart of non-subsampled Shearlet transform

1) 多尺度分解:利用非下采样金字塔(non-sampling pyramid, NSP)滤波器组进行尺度滤波,经过$j$层金字塔分解,可以得到$j$幅高频细节图像和低频近似图像,各子带分解图像大小与原图相同;

2) 方向剪切分解:运用剪切滤波器对高频图层图像进行方向分解,经过$M$级方向分解可得到$n = {2^M} + 2$个方向子带图像。

2.2 梯度绝对值能量算子

利用边缘梯度函数数值来判定开展图像的清晰度是典型的评价方法,常用的评价函数包括Tenengrad函数、方差(Variance)函数、能量梯度函数、Brenner函数、拉普拉斯函数、Robert函数等。理想的清晰度评价算子应该能够较好地反映图像的边缘特征,表征不同图像的清晰度差异,同时也应具有单峰性、灵敏性、抗噪性等性质。尤其是当图像中存在噪声时,会使得评价方法存在局部的波动,出现次峰,影响评价精度。在NSST变换域,本文利用梯度绝对值能量算子来构造清晰度评价算法。频域子带系数梯度绝对值能量和算子(sum of absolute gradient, SAG,用ESAG表示)定义为

${E_{{\rm{SAG}}}}(n, k) = \sum\limits_i {\sum\limits_j {A{G_{n, k}}(i, j)} } , \;\;A{G_{n, k}}(i, j) > T, $ (4)
$A{G_{n, k}}(i, j) = |{K_{n, k}}(i + 1, j) - {K_{n, k}}(i, j)| \\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, + |{K_{n, k}}(i, j + 1) - {K_{n, k}}(i, j)|, $ (5)

其中:${K_{n, k}}(i, j)$为待评价图像经过NSST分解后的子带系数,T为抗噪阈值。

SAG算子采用绝对值计算代替了能量梯度函数中的平方算子,相比于能量梯度函数中采用梯度平方和来作为评价数值,SAG算子在数学运算中将乘法运算操作转换为符号判断操作,在表征图像边缘清晰度的同时,降低了计算的复杂度,提高了运算效率。此外,引入了抗噪阈值,在计算评价函数值时只计算大于阈值的频域子带数值,增强了算法的抗噪性能。

2.3 显微图像清晰度评价

本文设计的结合多尺度分解和梯度绝对值算子的显微图像清晰度评价框架流程图如图 2所示。利用NSST对输入的待评价显微图像进行多尺度分解,得到1幅低频子带图像和若干幅不同尺度下不同方向的高频子带图像。计算各子带图像的梯度绝对值能量和算子,并通过设置抗噪阈值来排除噪声对评价结果的干扰,进而得到加权的高频子带能量和算子以及低频子带能量和算子,最终通过两者的比值得到清晰度评价算子。

图 2 结合多尺度分解和梯度绝对值算子的显微图像评价方法 Fig. 2 Micro-image definition evaluation method using multi-scale decomposition and gradient absolute value

NSST分解后的低频和高频子带信息具有不同的物理意义,高频信息对应图像的细节信息,低频信息对应图像的概貌信息,多尺度分解后的频域子带能够为图像清晰度评价提供依据。当显微图像出现离焦模糊时,低频子带系数能量变化不明显,对于图像清晰度的改变不是很敏感;而高频子带系数能量对图像清晰度变化非常敏感,清晰度下降时高频子带系数能量也发生明显的下降。利用两者的变化特点,可以构建有效的清晰度评价值。

设待评价的显微图像为I,所提出的图像清晰度评价算法处理步骤为

1) 对I进行非下采样剪切波多尺度分解,得到基础图层低频子带系数图像和若干不同方向的细节图层高频子带系数图像,非下采样多尺度分解层数为$N$,对应每层的分解方向个数为${D_n}(n = 1, 2, .., N)$

2) 对于步骤1)得到的不同尺度、不同方向的子带频域系数,定义$K_N^{\rm{L}}(i, j)$为在坐标位置(i, j)的低通子带系数,$K_N^{\rm{H}}(i, j)$n尺度k方向上在坐标位置(i, j)的各高频带通子带系数,其中$n = 1, 2, .., N$表示对应的分解尺度层数,k为剪切波方向子带的数目。

3) 计算得到基础图层的低频子带系数图像的能量和算子${R^{\rm{L}}}$

${R^{\rm{L}}} = E_{{\rm{SAG}}}^{\rm{L}}(N)。$ (6)

计算得到细节图层的高频子带图像的加权总能量和算子${R^{\rm{H}}}$

${R^{\rm{H}}} = \sum\nolimits_{n = 1}^N {{\gamma _n}R_n^{\rm{H}}} ,$ (7)

其中:$\sum \gamma = 1$, ${\gamma _n}$为各层子带的加权系数,表示不同尺度高频子带系数评价值的权重大小,通常尺度层数越小,设置权重越大,即${\gamma _1} > {\gamma _2} > ... > {\gamma _N}$$R_n^{\rm{H}}$为各细节尺度能量,通过计算该尺度下不同方向的高频子带图像能量平均值获得,表示为

$R_n^{\rm{H}} = \frac{{\sum {R_{n, k}^{\rm{H}}} }}{{{D_n}}},$ (8)

其中:Dn为对应尺度的分解方向个数,$R_{n, k}^{\rm{H}}$为由对应尺度的高频带通方向子带系数计算得到的各方向高频子带图像能量:

$R_{n, k}^{\rm{H}} = [{E_{{\rm{SAG}}}}(n, 1), {E_{{\rm{SAG}}}}(n, 2), \ldots , {E_{{\rm{SAG}}}}(n, k)]。$ (9)

4) 结合低频和高频子带能量和算子,最终的显微图像清晰度评价数值$C$由高频子带图像加权的总能量和算子除以低频子带图像的能量和算子计算得到,表示为

$C = \frac{{{R^{\rm{H}}}}}{{{R^{\rm{L}}}}}。$ (10)

当显微图像出现离焦模糊时,低频子带系数能量变化不明显,其对于图像清晰度的改变不是很敏感;而高频子带系数能量对图像清晰度变化非常敏感,图像清晰度下降时,高频子带系数能量也发生明显的下降。所提出的清晰度评值算子,评价值越大,表明图像清晰度越好,图像质量越好。

3 实验与分析 3.1 算法参数设置与说明

本文实验中,NSST分解层数$N = 3$,即将图像分解为1个低频子带和3个高频子带。对于3个高频子带,每个高频子带的分解方向数目分别为8,8和16,不同的方向代表了剪切波支撑基的不同剪切方向,可以解析出图像不同方向的细节信息。每层方向滤波器的大小分别为32×32,32×32和16×16。

在各子带系数图像梯度绝对值能量和的计算中,引入了噪声阈值,只计算大于噪声阈值的梯度能量,显著增强了算法抗噪性表现,实验中T取值范围为[1.5, 2],可根据图像内容和敏感度要求调节设置。低频子带图像的能量和由低频子带图像计算得到,各尺度高频子带图像的能量和是通过各尺度下不同方向的高频子带图像能量平均值获得。总的高频子带图像能量和是由各尺度高频子带能量加权求和得到,加权系数设置为${\gamma _1} = 0.7$${\gamma _2} = 0.2$${\gamma _3} = 0.1$。算法参数设置为经验值,可以根据具体应用场景和图像内容适当调整。

3.2 对比方法

本文引入多种清晰度评价方法,开展算法效果对比分析实验,对比方法包括Tenengrad函数[7]、能量梯度(energy of gradient, EOG)函数[10],拉普拉斯算子和(LS)[11],灰度平均梯度(grayscale mean gradient, GMG)算法[9]

Tenengrad函数用Sobel算子来估计水平方向和垂直方向梯度值,利用平方运算来表示图像边缘梯度大小。对于大小为M×N的待评价图像W,Tenengrad函数定义梯度平方和,表示为

$ T = \frac{1}{{(M - 2)(N - 2)}} \cdot {\sum\limits_i {\sum\limits_j {\left[ {S\left( {i, j} \right)} \right]} } ^2}, $ (11)

式中:$S(i, j)$表示图像$W$中像素$(i, j)$与Sobel算子的卷积,表示为

$S(i, j) = \sqrt {G_x^2(i, j) + G_y^2(i, j)}。$ (12)

而数字图像中,${G_x}(i, j)$${G_y}(i, j)$进而可以用下式表示:

$\begin{gathered} {G_x}(i, j) = W(i - 1, j + 1) + 2W(i, j + 1) \\ + W(i + 1, j + 1)\; - W(i - 1, j - 1) \\ - 2W(i, j - 1) - W(i + 1, j - 1), \end{gathered} $ (13)
$\begin{gathered} {G_y}(i, j) = W(i - 1, j - 1) + 2W(i - 1, j) \\ + W(i - 1, j + 1) - W(i + 1, j - 1) \\ - 2W(i + 1, j) - W(i + 1, j + 1), \end{gathered} $ (14)

其中:${G_x}(i, j)$${G_y}(i, j)$分别为图像W在像素$(i, j)$处的水平方向和垂直方向梯度评估数值。

能量梯度函数(EOG,用EEOG表示)利用相邻点的差分来计算一个点的梯度值,其评价算子表示为

${E_{{\rm{EOG}}}} = \frac{1}{{(M - 1)(N - 1)}}\sum\limits_i {\sum\limits_j {\{ {{[W(i + {\rm{1}}, j) - W(i, j)]}^{\rm{2}}}} } \\ + {[W(i, j + 1) - W(i, j)]^{\rm{2}}}\} 。$ (15)

拉普拉斯算子和(LS)比较了3×3像素邻域内的像素偏差程度,表示为

${A_{{\rm{LS}}}} = \frac{1}{{(M - 2)(N - 2)}} \\ \cdot \sum\limits_{i = 2}^{M - 1} {\sum\limits_{j = 2}^{N - 1} {\left| {8W(i, j) - W(i - 1, j - 1) - W(i - 1, j)} \right.} } \\ - W(i - 1, j + 1) - W(i, j - 1) - W(i, j + 1) \\ \left. { - W(i + 1, j - 1) - W(i + 1, j) - W(i + 1, j + 1)} \right|。$ (16)

LS评价值越大,说明图像信息越丰富,边缘越锐利,图像质量越好。

灰度平均梯度(GMG,用GGMG表示)评价指标反映了图像的细节和纹理变化特征,利用图像灰度梯度的平均值来表征图像的清晰程度,定义为

${G_{{\rm{GMG}}}} = \frac{1}{{(M - 1)(N - 1)}} \\ \cdot \sum\limits_{i = 1}^{M - 1} {\sum\limits_{j = 1}^{N - 1} {\sqrt {\frac{{{{[W(i + 1, j) - W(i, j)]}^2} + [{{(W(i, j + 1) - W(i, j)]}^2}}}{2}} } } 。$ (17)

GMG评价数值越大,表明纹理越丰富,图像越清晰。

此外,为了更清晰地对比本文方法和空域梯度绝对值能量函数(absolute gradient, AG,用GAG表示)的性能表现,将空域AG函数也作为对比算法,其定义为

${G_{{\rm{AG}}}} = \frac{1}{{(M - 1)(N - 1)}} \\ \cdot \sum\limits_i {\sum\limits_j {(\left| {W(i + 1, j) - W(i, j)} \right| + \left| {W(i, j + 1) - W(i, j)} \right|)} } 。$ (18)

AG函数评价值越大,表明图像边缘梯度越大,图像越清晰。

3.3 仿真实验结果

本节选取一幅清晰显微图在Matlab中加入不同程度的高斯模糊,仿真得到不同清晰度的显微图像序列,图像大小为800$ \times $800。所采用的清晰显微图像如图 3(a)所示,图 3(b)~ 3(f)分别为在图 3(a)清晰图的基础上加入不同的高斯模糊得到的一系列模糊显微图像,高斯模糊核半径为15,标准差$\sigma $分别为1,2,3,4,5。高斯模糊核标准差越大,得到的图像越模糊,清晰度越差。

图 3 不同模糊程度的显微图像仿真结果。(a)清晰图;(b)高斯模糊图σ=1; (c)高斯模糊图σ=2; (d)高斯模糊图σ=3; (e)高斯模糊图σ=4; (f)高斯模糊图σ=5 Fig. 3 Simulation results for microscopic images of different image blur.(a) Sharp image; (b) Gaussian blurred image σ=1; (c) Gaussian blurred image σ=2; (d) Gaussian blurred image σ=3; (e) Gaussian blurred image σ=4; (f) Gaussian blurred image σ=5

利用3.2小节中的几种对比评价算法和本文所提出的算法对图 3中不同模糊程度的序列显微图像进行清晰度评价,评价对比结果如表 1所示,所有算法都取4位有效数值进行对比。可以看到,随着图像清晰度的下降,图像越来越模糊,几种对比算法对应的评价数值也都随之变小,体现了图像清晰度变化的趋势,具有较好的单调性能。进一步分析,本文所提出的评价方法,对于不同清晰度显微图像评价值数值相对差异更加明显,区分度更高。为了能够更为客观地表述不同算法评价数值的区分灵敏程度,引入了灵敏度参数进行说明。

表 1 几种对比方法对于仿真显微图像序列的评价结果 Table 1 Evaluation results for simulated microscopic image sequence of several compared methods
模糊核标准差 Tenengrad EOG LS GMG AG Proposed
σ=0 3.749x103 87.43 9.257 2.111 3.771 1.749x10-1
σ=1 2.484x103 42.48 3.536 1.472 2.585 1.065x10-1
σ=2 1.598x103 25.67 1.753 1.144 2.228 5.082x10-2
σ=3 1.189x103 18.82 1.163 0.9793 2.107 1.823x10-2
σ=4 1.014x103 15.97 0.9316 0.9022 2.059 4.705x10-3
σ=5 9.356x102 14.72 0.8342 0.8661 2.037 1.148x10-3

理想的图像清晰度评价函数应该具有较高的灵敏度,特别是在清晰图像评价峰值附近有较为尖锐的评价值梯度变化[15]。为了衡量评价数值变化的灵敏程度,将序列图像评价的最大值与最小值之差与评价最大值的比值,作为评价方法灵敏度参数评价因子:

$S = \frac{{{V_{\max }} - {V_{\min }}}}{{{V_{\max }}}}, $ (19)

其中:${V_{\max }}$${V_{\min }}$分别为序列图像评价数值的最大值和最小值,$S$为灵敏度评价参数,$S$是介于0到1之间的数值,评价参数越大表示评价函数灵敏度越高。对于图 3序列图像的评价结果,各种对比方法的灵敏度参数如表 2所示。从表 2数据可以看到,所提出的评价方法具有较高的灵敏度,其对于不同清晰度图像的评价值区别较好。

表 2 几种清晰度评价方法灵敏度参数比较 Table 2 Comparisons of sensitivity parameter for several definition evaluation methods
评价方法 Tenengrad EOG LS GMG AG Proposed
灵敏度数值 0.750 0.832 0.910 0.590 0.460 0.994

为了进一步说明算法的抗噪特性的性能表现,对于图 3(a)中的清晰显微图像,加入不同程度的高斯模糊,高斯标准差$\sigma $变化范围为0~6,得到25张不同清晰度的序列图像。同时,在该序列图像中都加入高斯白噪声,高斯噪声均值为0,标准差为0.025,按照离焦模糊-聚焦清晰-离焦模糊的顺序排列,模拟显微图像自动对焦的过程。以其中几个位置的图像为例,含噪声的不同清晰度的显微图像序列示意图如图 4所示。其中,中间一幅为清晰图像,越往两侧图像模糊程度越大。

图 4 含噪声的不同清晰度的显微图像序列示意图 Fig. 4 The schematic diagram of noisy microscopic image sequences with different sharpness

对于上述带噪声的显微图像序列,利用几种对比的清晰度评价算法进行图像清晰度指标评价,得到归一化的评价结果曲线,如图 5所示。从图 5分析看到,归一化评价曲线基本反映了图像序列的清晰度变化趋势,但Tenengrad函数、LS评价算法的抗噪能力较差,评价曲线出现了明显的波动,形成了明显的次峰,影响了算法的准确性。虽然次峰的强度只有主峰能量的10%左右,但相对于次峰位置临近位置的评价值,次峰的扰动影响较大,造成了评价结果的无偏性和稳定性出现严重偏差。同时,相比于其他几种评价算法,LS算法的评价曲线半高宽度较窄,意味着评价函数的有效覆盖较小,当模糊程度增大时,算法评价值区分度较差。而本文所提出算法抗噪能力较强,基本不受噪声的干扰,表现出了较好的稳定性和鲁棒性,算法评价曲线变化相对缓慢,评价值覆盖范围更广。

图 5 噪声显微图像序列归一化评价结果曲线 Fig. 5 Normalized evaluation value curves for noisy microscopic image sequencesNormalized evaluation value
3.4 实拍实验结果

为了进一步验证评价算法对于实际显微镜拍摄图像的应用效果,使用了一款Olympus CX 23系列显微镜进行实拍实验。所用显微镜电子目镜采集设备的图像像素大小为2592×1944,约520万像素,其中,传感器尺寸为5.70 mm×4.28 mm,对角线长度为1/2.5英寸(1英寸=2.54 cm,约7.13 mm),像元尺寸为2.2 m。针对两组显微样品标本,在10倍放大倍数的情况下,分别手动调节采集了11幅和10幅物距从小到大的测试样本图像,代表从离焦到清晰再到离焦的对焦过程,两组样本图像分别为植物细胞有丝分裂图和洋葱鳞片表皮细胞图,分别如图 6(a)图 6(b)所示。其中,选了像素大小为800×800的图像块内容作为清晰度算法评价窗口,如图 6中红色方框区域所标记区域。实际应用中,可以根据需求选取任意位置的合适大小和图像内容的图像块,作为清晰度评价算法的输入。

图 6 实拍显微样本图像。(a)植物细胞有丝分裂图;(b)洋葱鳞片表皮细胞图 Fig. 6 Actual micro-image samples.(a) Plant cell mitosis image; (b) Onion scale epidermal cell image

利用几种对比评价算法对图 6实拍显微图像序列进行清晰度评价,评价结果数据如表 3表 4所示,由于各评价方法结果绝对大小有所区别,均采用4位有效数字的形式进行表示。对于两组图像,由于图像细节内容有所不同,评价数值相对大小有所差异。而对于同一组序列图像,各评价方法的数值结果能够体现图像由离焦到清晰再到离焦的过程,分别在第6幅图和第5幅图时达到了评价结果的最大值,这和人眼主观清晰度判断是一致。为了能够更好地展示不同算法的评价结果,绘制得到了归一化的评价结果曲线图,如图 7所示。从图 7曲线可以看到,本文方法得到的评价曲线对于两组实拍图像序列评价表现都较为稳定,峰值附近曲线变化趋势更为明显,鲁棒性和灵敏度更好。而在图 6(a)植物细胞有丝分裂序列图评价中,AG算法出现峰值错误,与人眼清晰度主观评价出现了偏差。

表 3 植物细胞有丝分裂图像评价结果 Table 3 Evaluation results of plant cell mitosis image
图像序列 Tenengrad EOG LS GMG AG Proposed
1 122.9 3.717 4.720 0.4354 2.410 2.147x10-3
2 235.4 5.925 5.227 0.5497 2.891 4.465x10-3
3 358.5 8.412 5.706 0.6549 3.386 1.368x10-2
4 658.1 14.56 6.719 0.8618 4.455 2.166x10-2
5 1246 27.40 8.824 1.1820 5.860 4.892x10-2
6 1493 33.30 10.04 1.3031 5.749 5.367x10-2
7 664.9 14.97 7.062 0.8736 4.066 1.530x10-2
8 393.3 9.239 5.968 0.6864 3.369 7.278x10-3
9 255.0 6.388 5.370 0.5708 2.916 5.278x10-3
10 183.1 4.926 5.026 0.5012 2.641 2.070x10-3
11 136.7 4.000 4.806 0.4516 2.443 2.061x10-3

表 4 洋葱表皮细胞图像评价结果 Table 4 Evaluation results of onion scale epidermal cell image
图像序列 Tenengrad EOG LS GMG AG Proposed
1 299.0 7.304 5.494 0.6103 4.862 2.438x10-3
2 477.9 11.00 6.198 0.7489 5.245 5.652x10-3
3 587.6 13.36 6.642 0.8255 5.627 1.137x10-2
4 767.2 17.37 7.347 0.9411 6.169 3.526x10-2
5 919.7 20.99 8.023 1.035 6.551 4.223x10-2
6 732.7 16.89 7.368 0.9279 6.001 3.263x10-2
7 539.7 12.60 6.604 0.8017 5.379 1.764x10-2
8 381.5 9.156 5.936 0.6832 4.816 6.008x10-3
9 300.5 7.470 5.588 0.6172 4.493 2.877x10-3
10 231.7 6.027 5.276 0.5544 4.178 1.896x10-3

图 7 实拍显微图像序列归一化评价结果曲线。(a)植物细胞有丝分裂图;(b)洋葱鳞片表皮细胞图 Fig. 7 Normalized evaluation value curves for actual micro-image samples.(a) Plant cell mitosis image; (b) Onion scale epidermal cell image

利用3.3小节所述的灵敏度评价参数,对两组实拍图像序列评价结果计算其灵敏度,结果数据如表 5所示。从表 5数据可看到,本文方法相比于其他几种对比方法,灵敏度参数数值最高,并且对于两组实拍序列都有较为稳定的表现,表明算法性能不依赖于图像内容。而几种对比评价算法对于洋葱鳞片表皮细胞图灵敏度出现了明显的下降,算法表现稳定度不够好。

表 5 实拍显微图像清晰度评价方法灵敏度参数比较 Table 5 Comparisons of sensitivity parameter for several definition evaluation methods of actual micro-images
评价方法 Tenengrad EOG LS GMG AG Proposed
植物细胞有丝分裂图 0.918 0.889 0.530 0.666 0.589 0.958
洋葱鳞片表皮细胞图 0.748 0.713 0.342 0.464 0.362 0.955
4 结论

本文提出了一种针对显微图像的清晰度评价算法,结合多尺度分解工具和频域梯度绝对值算子的设计,实现了有效的显微图像清晰度质量评价。利用NSST的多方向、非下采样、平移不变的特性,对输入的显微图像进行多尺度分解,得到一幅低频子带图像和若干幅高频子带图像。计算各子带图像的梯度绝对值算子和,并设置了抗噪阈值函数,在计算评价函数值时只计算大于阈值的频域子带数值,增强了算法的抗噪性能。利用图像清晰度变化对于低频和高频子带系数影响的差异,将高低频梯度绝对值算子的比值作为最终的显微图像清晰度评价数值。开展了仿真论证实验和实拍论证实验,获取了清晰度不同的显微图像序列,利用几种对比价算法进行清晰度指标计算。实验结果表明,所提出的清晰度评价算法具有较好的单调性和抗噪特性,和几种经典的评价算法相比,本文方法得到的评价结果在灵敏度、稳定性和鲁棒性方面表现更为优异,具有很好的实际应用价值,对数字显微镜自动对焦和显微镜成像质量评价提供了有效的支持。

在下一步的研究工作中,将重点研究所提出的清晰度评价算法与显微成像自动对焦应用的结合,完善算法的性能表现。

参考文献
[1]
Ilhan H A, Dogar M, Ozcan M. Digital holographic microscopy and focusing methods based on image sharpness[J]. Journal of Microscopy, 2014, 255(3): 138-149. [Crossref]
[2]
Rudnaya M E, Mattheij R M M, Maubach J M L. Evaluating sharpness functions for automated scanning electron microscopy[J]. Journal of Microscopy, 2010, 240(1): 38-49. [Crossref]
[3]
Zheng X, Ai L F, Liu K, et al. Auto-focusing function for microscopic images based on global and local gray-scale variation[J]. Laser & Optoelectronics Progress, 2017, 54(8): 081801.
郑馨, 艾列富, 刘奎, 等. 结合全局和局部灰度变化的显微图像自动聚焦函数[J]. 激光与光电子学进展, 2017, 54(8): 081801 [Crossref]
[4]
Zhao J F, Mao L, Liu C, et al. Microscopy imaging definition criterion using visual attention mechanism and edge spreading evaluation[J]. Acta Photonica Sinica, 2015, 44(7): 0711002.
赵巨峰, 毛磊, 刘承, 等. 视觉注意机制与边缘展宽衡量相结合的显微成像清晰度评价[J]. 光子学报, 2015, 44(7): 0711002 [Crossref]
[5]
Zhou H K, Ge P S, Feng H L. Micro-image definition evaluation algorithm based on non-subsampled Contourlet transform[J]. Journal of Hefei University of Technology (Natural Science), 2014, 37(10): 1204-1209.
周厚奎, 葛品森, 冯海林. 基于非下采样Contourlet变换的显微图像清晰度评价算法[J]. 合肥工业大学学报(自然科学版), 2014, 37(10): 1204-1209 [Crossref]
[6]
Jiang M S, Zhang N N, Zhang X D, et al. Applications of hybrid search strategy in microscope autofocus[J]. Opto-Electronic Engineering, 2017, 44(7): 685-694.
江旻珊, 张楠楠, 张学典, 等. 混合搜索法在显微镜自动对焦中的应用[J]. 光电工程, 2017, 44(7): 685-694 [Crossref]
[7]
Huang D T. Study on auto-focusing method using image technology[D]. Changchun: Graduate School of Chinese Academy of Sciences (Changchun Institute of Optical Precision Machinery and Physics), 2013.
黄德天.基于图像技术的自动调焦方法研究[D].长春: 中国科学院研究生院(长春光学精密机械与物理研究所), 2013. [Crossref]
[8]
Wang P J, Chen D J, Wu M L, et al. Implementation of a multiple automatic focusing algorithm for video measuring system[J]. China Mechanical Engineering, 2007, 18(21): 2555-2560.
王平江, 陈德军, 巫孟良, 等. 一种复合的自动对焦方法在影像测量仪中的应用[J]. 中国机械工程, 2007, 18(21): 2555-2560 [Crossref]
[9]
Jin H Y, Wang Y Y. A fusion method for visible and infrared images based on contrast pyramid with teaching learning based optimization[J]. Infrared Physics & Technology, 2014, 64: 134-142. [Crossref]
[10]
Guo J B, Feng H J, Wang L, et al. Design of focusing window based on energy function of gradient[J]. Infrared Technology, 2016, 38(3): 197-202.
郭敬滨, 冯华杰, 王龙, 等. 基于梯度能量函数的调焦窗口构建方法[J]. 红外技术, 2016, 38(3): 197-202 [Crossref]
[11]
Yao S, Cao D H, Wu Y B. Performance of autofocus judgements in different quantity of information series images[J]. Opto-Electronic Engineering, 2006, 33(5): 81-84, 90.
姚松, 曹丹华, 吴裕斌. 图像信息量的变化对自动对焦评价函数的影响[J]. 光电工程, 2006, 33(5): 81-84, 90 [Crossref]
[12]
Zhang Y T, Ji S P, Wang Q F, et al. Definition evaluation algorithm based on regional contrast[J]. Journal of Applied Optics, 2012, 33(2): 293-299.
张亚涛, 吉书鹏, 王强锋, 等. 基于区域对比度的图像清晰度评价算法[J]. 应用光学, 2012, 33(2): 293-299 [Crossref]
[13]
Jin W Z, Leng Q J, Zhang Z, et al. Contourlet transform based multiscale image quality assessment metric[J]. Journal of Wuhan University (Natural Science Edition), 2015, 61(2): 192-196.
金伟正, 冷秋君, 张卓, 等. 基于Contourlet变换的多尺度图像质量评价[J]. 武汉大学学报(理学版), 2015, 61(2): 192-196 [Crossref]
[14]
Wu W, Qiu Z M, Zhao M, et al. Visible and infrared image fusion using NSST and deep Boltzmann machine[J]. Optik, 2018, 157: 334-342. [Crossref]
[15]
Wang Y R, Feng H J, Xu Z H, et al. An adjustable coverage range autofocus evaluation function using gradient operator with variable frequency[J]. Infrared and Laser Engineering, 2016, 45(10): 1028001.
王烨茹, 冯华君, 徐之海, 等. 一种覆盖范围可调的变频梯度自动对焦评价函数[J]. 红外与激光工程, 2016, 45(10): 1028001 [Crossref]