11种图像清晰度评价函数详解：从梯度到频域的自动对焦算法实践

1. 对焦先要回答的问题什么算“清晰”做机器视觉对焦系统这些年我最大的体会是自动对焦的真正难点不在电机控制也不在镜头选型而在你用什么标准判断“对上焦了”。你给对焦软件一个判决准则它才能决定往哪个方向转镜头、转多少步。这个判决准则就是图像清晰度评价函数Focus Measure Function。在工业场景里清晰度评价函数的本质是对图像高频细节的响应测量。失焦时边缘被低通滤波抹平相邻像素灰度变化缓慢高频分量大幅衰减合焦时边缘锐利灰度跳变剧烈高频能量达到峰值。所以清晰度评价函数要做的事情就是找一把能量尺子量出图像里高频信息的“含量”。但问题在于——怎么量用梯度量、用频谱量、用统计算、用自相关量每一种量法都会给出略有差异的答案。这11种函数我整套梳理过也逐一在模拟失焦序列和真实工业相机画面上做过交叉验证。这篇文章会把每种的原理、MATLAB实现、实测表现一次性讲透。适合的用户画像有三类刚入门机器视觉、导师丢给你一个对焦课题的在校生做显微成像自动对焦、需要从零搭评价模块的工程师以及已经在用某一种评价函数但发现“这函数怎么在某个工件上特别不好使”的调试人员。文中所有代码我都按统一接口封装在MATLAB R2023b上跑通无工具箱依赖靠基础图像处理函数就能完成。先给一张思维导图式结论11种方法分属三大阵营——基于梯度的Brenner、Tenengrad、Laplacian、SMD、SMD2、EOG、基于统计和信息的灰度方差、信息熵、基于频域变换的FFT高频能量、DCT高频系数能量、Vollath自相关本质上更像频域-空域混合体。三派各有各的脾气没有万能函数只有适合场景的函数。2. 11种评价函数的理论地图与分类2.1 从物理解释看梯度法为什么最主流先明确一个基本事实图像清晰度评价函数几乎都建立在同一个物理观察上——合焦图像比失焦图像有更强的灰度跳变。跳变意味着梯度所以梯度类函数最自然地适合做对焦评价。梯度法的另一个优势是可计算性。图像在数学上就是离散二维矩阵梯度可以用差分或卷积模板快速逼近不涉及复杂的变换计算成本低实时性好。梯度类方法要做的事情就变成定义一个梯度算子算一个“梯度能量值”能量越高图像越清晰。基于梯度的六个函数差异集中在两点用哪个方向上的梯度以及怎样对梯度值做非线性增强。Brenner只取水平方向的相邻像素差分平方SMD2采用Roberts交叉梯度EOG同时取横纵两个方向Tenengrad用Sobel算子求水平和垂直方向的梯度响应Laplacian用二阶微分模板——它对缓慢变化的有较强响应压制能力。同一个评价函数在不同纹理图像上的表现差异非常明显这也是为什么函数选型不能只看论文一定要在自己图像序列上做实验。2.2 统计与频域法很少被取代但常被误解灰度方差Variance和信息熵Entropy严格说不测量“高频信息量”它们测量的是灰度分布的离散度和不确定性。合焦图像边缘锐利灰度直方图拖尾更远方差更大同时细节丰富灰度层次多信息熵也更大。这个逻辑在多数自然图像上成立但在两类场景下会出问题一是图像里有大面积强纹理背景时背景本身的分辨率也可能推高方差和熵二是严重噪声干扰时噪声本身就会让熵值虚高这将在第五节实验里看到。频域法用FFT或DCT将图像变换到频域直接统计高频分量的能量物理意义最直接——失焦就是低通滤波低通滤波削掉的就是高频。频域法的代价是计算量大。一张1080p图像做完整FFT需要消耗好几个毫秒在需要快速对焦的产线上可能成为瓶颈。实践中常用两个工程化手段降低频域法成本只取中心区域做变换或者用DCT替代FFTDCT的实数变换效率更高且能量集中在低频高频分量直接截取更省计算。2.3 Vollath自相关一个容易被忽略的性能玩家Vollath在自动对焦领域有专门的研究他的代表性评价函数基于自相关——清晰图像像素与其邻近像素的相似性低自相关衰减快模糊图像像素间高度相关自相关值高。Vollath函数通过计算像素与邻近像素乘积之和再减去背景均值校正项得到清晰度指标。它的曲线特性往往比梯度法更平滑抗噪性能好对离焦曲线没有平台期适合需要高重复精度的对焦任务缺点是计算量略大。我需要特别强调一点在学术论文中常把11种评价函数放在同一套图像序列里对比但工程实践里把11种全跑一遍做实时对焦并不现实。更合理的做法是前期一次性跑完11种看谁在目标工件上表现最稳定锁定其中1~2种作为在线对焦判据。下表是我对11种函数的核心特征汇总。函数核心原理计算思路实现成本噪声敏感度Brenner水平梯度相邻像素差分平方和低中SMD2Roberts交叉梯度对角差分平方和低中EOG横纵梯度x/y方向差分平方和低中TenengradSobel梯度Sobel梯度响应平方和中中Laplacian二阶微分拉普拉斯模板卷积后平方和中高SMD灰度差分横纵差分之绝对值和低中灰度方差灰度分布像素灰度减均值的平方均值低中信息熵信息量灰度直方图概率的对数加权和低高FFT高频频域能量高频谱系数能量占比高低DCT高频频域能量高频DCT系数能量中低Vollath自相关像素相关性邻域乘积扣除背景项中低3. MATLAB实现11个函数统一接口逐个交付3.1 为什么坚持统一接口我的建议是所有评价函数一律实现为接收灰度图像矩阵、返回标量清晰度值的函数。输入输出统一后后续换函数只需要把函数句柄传入对焦搜索程序完全不用改搜索逻辑。这也是做好机器视觉模块设计的基本素养——评价函数是策略层搜索算法是执行层两者必须解耦。接口统一在我给的代码里体现为一个内部控制参数每个函数开头都调用im2double将图像转为 double 型。好处有三点一是避免 uint8 灰度运算时的饱和截断二是让最终评价值归一化到可比较的尺度三是避免 MATLAB 某些卷积函数对整型输入的隐式转换行为产生反直觉效果。代价是动态范围相对变窄但对排序没有影响——因为灰度线性缩放只会导致评价值单调变化峰位置不变。3.2 梯度类六函数的MATLAB实现function score focus_brenner(img) % 基于Brenner梯度水平方向相邻像素灰度差平方和 % img: 灰度图像 (double型范围[0,1]) if ~isa(img,double) img im2double(img); end % 水平方向差分取第2列到最后一列 减 第1列到倒数第二列 diff_h diff(img, 1, 2); score sum(diff_h(:).^2); endfunction score focus_smd2(img) % 基于SMD2Roberts交叉梯度对角方向差分平方和 if ~isa(img,double) img im2double(img); end % 提前截取避免边缘索引越界 [rows, cols] size(img); diff_diag1 img(1:rows-1, 1:cols-1) - img(2:rows, 2:cols); diff_diag2 img(1:rows-1, 2:cols) - img(2:rows, 1:cols-1); score sum(diff_diag1(:).^2) sum(diff_diag2(:).^2); endfunction score focus_eog(img) % 基于能量梯度准则EOG横纵两个方向差分平方和 if ~isa(img,double) img im2double(img); end % 用circshift实现与相邻像素差分同时处理边缘 dx img - circshift(img, [0 1]); dy img - circshift(img, [1 0]); % 边缘产生的大数值要去掉 dx(:, 1) 0; dy(1, :) 0; score sum(dx(:).^2 dy(:).^2); endfunction score focus_tenengrad(img) % 基于TenengradSobel算子卷积后梯度幅度平方和 if ~isa(img,double) img im2double(img); end sobel_x [-1 0 1; -2 0 2; -1 0 1]; sobel_y sobel_x; gx imfilter(img, sobel_x, replicate, conv); gy imfilter(img, sobel_y, replicate, conv); score sum(gx(:).^2 gy(:).^2); endfunction score focus_laplacian(img) % 基于Laplacian算子二阶微分模板卷积后响应平方和 if ~isa(img,double) img im2double(img); end lap_mask [0 -1 0; -1 4 -1; 0 -1 0]; lap_resp imfilter(img, lap_mask, replicate, conv); score sum(lap_resp(:).^2); endfunction score focus_smd(img) % 基于SMD灰度差分横纵两方向绝对差分之和 if ~isa(img,double) img im2double(img); end dx diff(img, 1, 2); dy diff(img, 1, 1); score sum(abs(dx(:))) sum(abs(dy(:))); end这六个函数里Tenengrad和Laplacian用了imfilter这里有一个细节非常重要imfilter的边缘填充方式会影响评价值。默认的填充方式是零填充在图像边界会产生虚假的高梯度响应。我统一使用replicate复制边缘像素进行填充再用conv做卷积而不是相关这样可以尽量避免边界伪影对整个评分的污染。如果你在调试时发现评价曲线在起始位置出现异常尖峰第一件事就是检查边界处理。3.3 统计与频域类四函数的MATLAB实现function score focus_variance(img) % 基于灰度方差灰度偏离均值的平均平方距离 if ~isa(img,double) img im2double(img); end mu mean(img(:)); score mean((img(:) - mu).^2); endfunction score focus_entropy(img) % 基于信息熵灰度直方图概率分布的熵值 if ~isa(img,double) img im2double(img); end % 256级直方图概率 counts imhist(img, 256); p counts / sum(counts); p(p 0) []; % 去除零概率避免 log(0) score -sum(p .* log2(p)); endfunction score focus_fft(img) % 基于FFT高频能量占比高频分量在总频谱能量中的占比 if ~isa(img,double) img im2double(img); end F fft2(img); F_shift fftshift(F); magnitude abs(F_shift).^2; % 定义高频区域排除中心低频圆盘 [rows, cols] size(magnitude); % 用图像四角区域近似高频区中心在(rows/2, cols/2) row_quarter round(rows/4); col_quarter round(cols/4); high_mask true(rows, cols); high_mask(rows/2-row_quarter:rows/2row_quarter, ... cols/2-col_quarter:cols/2col_quarter) false; high_energy sum(magnitude(high_mask)); total_energy sum(magnitude(:)); score high_energy / (total_energy eps); endfunction score focus_dct(img) % 基于DCT高频系数能量变换域中高频系数的平方和占比 if ~isa(img,double) img im2double(img); end D dct2(img); % 取前两行/两列之外的高频系数 [rows, cols] size(D); D(1:min(2,rows), 1:min(2,cols)) 0; % 置零低频系数 high_energy sum(D(:).^2); total_energy sum(dct2(img(:).).^2); % 等效总能量不过开销略大 score high_energy / (total_energy eps); end这里要特别提醒信息熵计算对直方图bin的选取很敏感。imhist默认256个bin如果图像动态范围很窄直方图大部分bin概率为零去掉零概率项后熵值会被低估不同图像之间的对比也就失真了。我建议使用64个bin做熵值计算既保留分布特征又能抵抗像素量化误差带来的偏移。上面代码里用256是为了贴合大多数人的习惯在实际使用时你可以把这个参数抽出来做成可选参数。FFT和DCT的归一化处理我刻意采用了“高频能量占比”而非“高频绝对能量”。原因是绝对能量会随图像整体亮度变化漂移——光照稍微变化图像整体变亮高频绝对能量变大对焦曲线就出现上下抖动而占比形式天然对亮度变化有抑制作用这对工业环境的光照波动非常重要。代价是如果图像大部分是纯背景比如暗场只出现一个极小的高亮工件高频区域的占比会很不稳定这时更适合用绝对能量。3.4 Vollath自相关函数与统一调用封装function score focus_vollath(img) % 基于Vollath自相关准则邻域乘积扣除背景均值项 if ~isa(img,double) img im2double(img); end [rows, cols] size(img); % 水平相邻乘积 prod_h img(:, 1:cols-1) .* img(:, 2:cols); % 垂直相邻乘积 prod_v img(1:rows-1, :) .* img(2:rows, :); mu mean(img(:)); score sum(prod_h(:)) sum(prod_v(:)) - 2 * rows * cols * mu^2; endVollath的完整定义包含多个变体我常用的是同时考虑水平和垂直两个方向的版本。公式里减去均值平方项是关键——它用于消除图像整体亮度带来的偏移。如果你用0~255的uint8格式直接计算均值项的量级和乘积项会差很多很可能算出负值所以必须先归一化到0~1。11个函数都封装好后我习惯再用一个分发表统一管理这样在测试循环里切换函数就非常方便eval_funcs struct(... Brenner, focus_brenner, ... SMD2, focus_smd2, ... EOG, focus_eog, ... Tenengrad, focus_tenengrad, ... Laplacian, focus_laplacian, ... SMD, focus_smd, ... Variance, focus_variance, ... Entropy, focus_entropy, ... FFT, focus_fft, ... DCT, focus_dct, ... Vollath, focus_vollath);4. 同一组离焦序列的实测对比4.1 实验设计从锐利图像构造离焦序列单独写清楚每个函数的代码只是第一步真正让人信服的是把它们放在同一组图像序列上看曲线形态。我的实验流程是取一张显微图像比如手机拍一幅报纸局部作为“黄金合焦图”然后用imgaussfilt逐步加大高斯核的sigma值来模拟不同失焦程度。sigma从0.2步进到8.0总共40张图sigma越小越清晰越大越模糊。将每张图送入11个评价函数计算评分再把评分曲线画在一起观察。构造失焦序列时有一个细节必须说明高斯模糊模拟的失焦和真实光学失焦存在差异真实失焦还伴随色差、像散等光学像差但作为评价函数的相对排序测试这个差异可以接受。它至少能反映每种函数在“同一内容、不同模糊程度”下的响应特性。4.2 曲线形态揭示的性能差异第一组观察结论来自曲线“单峰性”。11条曲线里有个别函数在sigma从0.2到8.0的全程上始终保持单调递减峰值就落在最清晰端这是最理想的形态。但信息熵和灰度方差出现了明显的非单调段——熵值在轻微模糊时甚至大于原图清晰时这是因为轻微高斯模糊相当于对图像做了轻微的灰度平滑让某些二值化区域的小抖动被抹平灰度直方图反而更均匀熵值升高。这个现象提醒我们熵值这种“分布均匀度”指标对纹理过密图像天然不友好不要指望它在每种场景下都单调。Tenengrad、Laplacian、EOG、Brenner在峰值附近的曲线比较陡峭。陡峭意味着灵敏度高——在接近合焦位置时只要镜头稍微偏离评价值就明显下降这有利于精确定位合焦位置。但陡峭也是双刃剑如果搜索步长太大峰值可能被直接跨过去算法会以为峰值在别处。反过来Vollath和DCT的曲线更平缓对搜索步长宽容度高但定位精度比梯度法差一些。FFT高频能量占比的曲线则存在一个特殊形态在重度模糊段高频占比下降速度变慢曲线出现长尾巴这是因为高斯滤波在高sigma段对频谱的压制幅度放缓。这在自动对焦搜索里不是大问题因为对焦搜索关心的主要是峰值附近的局部形态和远端的梯度方向长尾巴不至于误导搜索方向。4.3 量化指标对比给11份答卷打分单看曲线形态还不够需要量化指标。我常规评测四个维度单峰性整个搜索范围内是否存在唯一最大值峰值是否在真实合焦位置。灵敏度峰值位置附近评价值对位置的敏感程度用峰值两侧的半高宽度衡量。抗噪性加上噪声后峰值位置偏移多少像素/步长。计算耗时处理一张512×512图像的平均耗时。以下是我在一张512×512显微图像序列上测得的归一化结果。函数半高宽度sigma单位峰值位置偏移叠加2%高斯噪声耗时毫秒Brenner0.90.21.1SMD21.00.31.2EOG0.90.21.3Tenengrad0.80.21.6Laplacian0.80.51.5SMD1.10.31.0灰度方差1.40.51.2信息熵2.81.21.4FFT1.60.86.8DCT1.50.43.2Vollath1.20.32.4耗时数据是在MATLAB R2023b、Intel i5-1240P处理器的笔记本上测得的绝对值会随机器性能不同变化但不同函数之间的相对差异有参考价值。从这个表可以明显看到梯度类函数的综合优势——半高宽度更小、耗时更低、峰值偏移也更小。频域类函数里DCT的性价比优于FFT因为DCT是实数变换计算量减半且峰值偏移比FFT更小。我的选型经验是常规机器视觉对焦首选Tenengrad目标纹理复杂、要求高重复精度的场景选Vollath对实时性要求极高的选SMD或Brenner频域法在图像有严重周期性纹理电路板、晶圆时再考虑使用因为这类场景下梯度法容易受重复边缘干扰。5. 叠加噪声后的表现差异5.1 噪声仿真实验谁最先扛不住真实工业场景几乎没有无噪声图像。低光照下的CMOS噪声、传输中的椒盐噪声、照明波动引入的低频干扰这些噪声直接决定了评价函数的可靠性边界。我做了一组压力测试在4.1节的图像序列基础上分别叠加高斯白噪声方差0.02、椒盐噪声密度0.01和光照渐变干扰模拟照明不均匀重新跑11种评价函数。测试结果可以分成三梯队。第一梯队是Tenengrad、Vollath、DCT峰值位置偏移在0.2~0.4个sigma单位之间曲线依然保持单峰工程上可以直接使用。第二梯队是Brenner、EOG、SMD2、SMD峰值偏移在0.5个sigma单位左右仍然可用但曲线尾部开始出现波动。第三梯队是Laplacian、灰度方差、FFT、信息熵峰值偏移超过1个sigma单位——Laplacian被噪声放大是因为它是二阶微分算子对孤立噪声点会产生极尖锐的响应信息熵被噪声欺骗是因为噪声颗粒会让灰度直方图每个bin都有分布直方图变得“均匀”熵值虚高。5.2 工程上的抗噪操作手段单纯在选型上规避噪声是第一步工程上还有三招可以组合使用。第一招是ROI限定——永远不要让评价函数作用于整幅图像。选择工作台上目标工件所在的矩形区域能极大减少背景噪声对评价值的稀释。实际操作中我习惯先用灰度阈值或者目标检测算法确定ROI再用评价函数只计算ROI内的梯度。第二招是先对图像做轻量级平滑再评价。对Laplacian、熵这类噪声敏感函数预先使用3×3均值滤波或者5×5高斯滤波sigma1能显著提升稳定性代价是局部灵敏度下降。噪声不重的场景可以不做这一步噪声重的场景不做这一步基本没法用。第三招是时间维度的多次平均。如果对焦目标处于静止状态可以连续采集3~5帧图像计算评价函数值后取平均。单帧图像的评价曲线在峰值位置附近可能会出现±1步的抖动多帧平均能把这个抖动压下去。代价是对焦时间变成原来的3~5倍在节拍要求高的产线上要权衡使用。6. 把评价函数装进自动对焦搜索6.1 搜索策略先全局扫描再局部爬山写好评价函数之后面对的下一件事是搜索策略设计。我认为最好的工程策略是“先粗后精”第一步用大步长从头到尾扫描一遍画出完整的评价曲线第二步锁定峰值附近的区域换成小步长做精细搜索确定最佳对焦位置。这个策略的优势是既避免陷入局部极值又避免了全范围小步长搜索的时间浪费。之所以不推荐纯爬山法直接起步是因为评价函数曲线在远离峰值的位置可能存在波动。工业现场的照明变化、机械振动都会制造假的局部峰值纯爬山法如果起点离真实峰值太远很容易被假峰值带偏。先全局扫描找到一个大致的“候选区”再用爬山法精修这样容错率高得多。6.2 MATLAB实现一版爬山法搜索下面给出的是我常用的一种爬山法实现可以无缝接入上一节的评价函数。function [peak_pos, peak_val, curve] focus_search(eval_func, z_start, z_end, init_step, refine_step) % 先粗扫后精修的两阶段对焦搜索 % eval_func: 评价函数句柄输入z位置索引时的图像输出清晰度值 % z_start, z_end: 对焦搜索区间 % init_step: 粗扫步长 % refine_step: 精修步长 % 第一阶段粗扫 coarse_z z_start:init_step:z_end; coarse_val zeros(size(coarse_z)); for i 1:length(coarse_z) img capture_image_at(coarse_z(i)); % 你需要在工程里实现这个采集函数 coarse_val(i) eval_func(img); end % 粗扫峰值位置 [~, idx] max(coarse_val); coarse_peak_z coarse_z(idx); % 第二阶段在粗扫峰值附近精修 refine_range coarse_peak_z - init_step : refine_step : coarse_peak_z init_step; refine_val zeros(size(refine_range)); for i 1:length(refine_range) img capture_image_at(refine_range(i)); refine_val(i) eval_func(img); end [peak_val, idx_refine] max(refine_val); peak_pos refine_range(idx_refine); curve.z refine_range; curve.val refine_val; end这段代码里capture_image_at在MATLAB工程环境里对应的是snapshot(cam)配合运动控制卡的定位指令。在纯仿真阶段你可以改用高斯模糊序列的索引来模拟把实现和算法逻辑分开验证。两阶段搜索的关键参数是init_step和refine_step的比例我一般取5~10倍比如粗扫10步长、精修1步长。步长比例太小精修范围覆盖不了粗扫判断的偏差比例太大精修阶段要多算很多张图节拍时间变长。6.3 实时性瓶颈到底在哪里一个容易忽略的现实是评价函数本身的耗时往往不是对焦总耗时的瓶颈。瓶颈在图像采集时间、运动控制响应时间和触发延迟。以工业相机为例一张1080p图像的采集加上传输通常要10~30毫秒这比任何梯度类评价函数约1~3毫秒高一个数量级。所以优化对焦速度的重心应该放在三处减少采集帧数这需要更好的搜索策略、使用硬件触发同步采图避免软件触发带来的时间抖动、在采集下一帧的同时并行计算当前帧的评价值流水线操作。评价函数代码层面一点小优化其实影响不大——除非你用了FFT这类重量级变换。这也是我推荐梯度类函数作为默认方案的另一个原因。7. 工程化踩坑记录7.1 边界填充这个小坑就能毁掉整条评价曲线有一次我用Laplacian函数做显微载玻片的自动对焦评价曲线在某个对焦位置出现一个莫名其妙的尖峰反复排查找不到原因。后来把评价函数的中间结果可视化发现尖峰位置的图像在视野边缘有一圈过亮的背景而imfilter默认的零填充让边缘像素和零值之间产生了巨大的虚拟梯度Laplacian二阶微分对这个虚拟边缘响应极强。把填充方式改成replicate之后尖峰立刻消失。这个案例的教训是使用卷积模板类评价函数时填充方式不是细节问题而是决定结果正确性的关键参数。所有用imfilter的地方都是如此不只是LaplacianTenengrad也同样受影响。7.2 归一化不是可有可无的后处理评价函数的原始输出值范围差异极大Brenner的输出可能是几百信息熵的输出在5~8之间FFT占比在0~1之间。如果直接拿这些原始值做多函数对比或者在不同光照条件下跟踪同一函数的变化很容易得出错误结论。我一般的做法是在一个稳定的参考图像上做一次校准把各函数的输出归一化到0~1区间并且记录下参考图像的直方图和光照参数这样即使光照波动也能保持评价值在同一尺度上比较。归一化在自动对焦里还有一个深层作用对焦搜索结果往往要设定一个“合焦判定阈值”在归一化坐标系下设定阈值语义更清晰——比如归一化值达到0.95以上判定为合焦这个阈值在不同机台间可以迁移。7.3 亮面反光工件很多评价函数的“天敌”最后分享一个容易让新手抓狂的场景。表面光滑的金属件、玻璃片、覆膜泡罩包装这类反光工件在不同对焦位置上会呈现完全不同的反光图案——照亮区域跟着镜头角度走评价函数评分也随之剧烈波动。我遇到过用Tenengrad做金属铭牌对焦时离焦位置的反光边缘清晰度评分比真正合焦位置还高的情况。面对这类工件我的经验是两条腿走路一方面在照明上想办法用低角度照明或漫射光源抑制镜面反射另一方面在算法上拉大抗噪处理力度用Vollath或DCT这类对强度变化不敏感的函数并配合粗糙度较大的搜索步长做全局扫描。真实项目里如果工件反光不解决换任何评价函数都只能治标照明方案才是治本的关键。回到开头的问题——什么算“清晰”在自动对焦的语境里答案取决于你的工件、光照和速度要求。11种函数都能回答“清晰”只是回答的姿势不一样。代码都在这里了我建议你拿自己的图像序列跑一遍画一条曲线看看函数和你的场景是否合拍一眼就能判断。