Halcon手写Otsu最大类间方差法：从原理到产线级实现

1. 为什么要在Halcon里手写Otsu1.1 从一次产线调试的尴尬说起前阵子帮一个做五金件表面检测的朋友调一条老产线相机是五百万像素的卷帘快门光源用的是环形白光工件表面有拉丝纹路还有反光。他们原来的方案是用固定阈值做二值化白天和夜班切换的时候车间顶灯亮度不一样阈值就得手动改操作员改烦了干脆把阈值调得很宽结果缺陷漏检率直接飙上去。我当时第一反应是上Halcon自带的binary_threshold算子选max_separability一行代码就能出结果省事。但朋友那边有个硬性约束他们的视觉软件跑在一台老工控机上Halcon版本是几年前的老版本binary_threshold的max_separability模式在某些图像上表现不稳定而且他们希望把阈值计算过程完全掌握在自己手里方便后续做自适应调整和日志记录。这就逼着我回到最原始的路子自己实现一遍Otsu最大类间方差法。说实话Otsu这个算法本身不复杂大学数字图像处理课上都讲过但真正在Halcon环境里用HDevelop或者C接口把它写稳、写快、写得能上产线还是有不少细节要抠。比如Halcon的灰度直方图算子返回的数据结构怎么用、循环里怎么避免频繁调用算子、阈值算出来之后怎么和Halcon的threshold算子衔接、多通道图像怎么处理、ROI怎么限定统计范围这些都是实际写的时候才会碰到的问题。这篇文章就是把这套东西完整梳理一遍。不管你是刚接触Halcon的新手还是用了几年但一直调库函数没自己写过底层算法的老手都能从里面找到能直接抄作业的代码和踩坑经验。核心关键词就一个Halcon实现Otsu最大类间方差法围绕它把原理、算子选型、代码实现、性能优化、常见问题全部讲透。1.2 Otsu到底解决了什么问题用大白话说Otsu就是帮你在灰度图像里自动找一个“最合适的分界线”把图像分成前景和背景两拨。这个“最合适”的标准是让两拨之间的差异尽可能大同时每拨内部的差异尽可能小。数学上就是用类间方差来衡量类间方差最大的那个灰度值就是我们要的阈值。它特别适合那种前景和背景灰度分布比较分明、但整体亮度会漂移的场景。比如上面说的五金件检测工件本体灰度大概在80到120之间波动背景传送带灰度在180到200之间中间有个明显的谷。Otsu就能自动找到谷底那个位置不用人天天去调。但如果图像里前景背景面积严重失衡比如缺陷只占整幅图的千分之一那Otsu算出来的阈值就会偏向大面积的那一类这时候就得配合ROI或者先做预处理。Halcon本身提供了binary_threshold算子里面Method参数选max_separability就是Otsu的实现。但为什么还要自己写三个原因一是老版本Halcon或者某些授权模式下这个算子行为不一致二是自己写可以灵活控制统计范围、可以加权重、可以做多级阈值扩展三是调试的时候能把中间变量打出来看心里有底。下面就从原理开始一步步把它落地。2. Otsu最大类间方差法的原理拆解2.1 类间方差的数学推导假设一幅灰度图像的灰度级是0到L-1总共L个级别Halcon里常见的8位图就是L256。我们用n_i表示灰度值为i的像素个数N是总像素数那么每个灰度级出现的概率就是p_i n_i / N现在假设我们选了一个阈值t把图像分成两类C0是灰度0到t的像素C1是灰度t1到L-1的像素。那么两类的出现概率分别是ω0 Σ(i0到t) p_i ω1 Σ(it1到L-1) p_i 1 - ω0两类的平均灰度分别是μ0 Σ(i0到t) i * p_i / ω0 μ1 Σ(it1到L-1) i * p_i / ω1整幅图像的平均灰度是μT Σ(i0到L-1) i * p_i ω0 * μ0 ω1 * μ1类间方差定义为σ²_B ω0 * (μ0 - μT)² ω1 * (μ1 - μT)²把μT代入化简可以得到一个更常用的形式σ²_B ω0 * ω1 * (μ0 - μ1)²这个形式计算起来更方便因为只需要算ω0、μ0、μ1三个量。Otsu算法就是遍历所有可能的t找到让σ²_B最大的那个t作为最终阈值。2.2 为什么类间方差最大就是最优这里稍微展开说一下直觉。σ²_B衡量的是两类均值之间的差距乘以ω0*ω1相当于给这个差距加了一个权重。当阈值取得太靠左C0里像素很少ω0很小乘积就小阈值太靠右ω1很小乘积也小。只有阈值落在两类分布中间的谷底附近ω0和ω1都不太小同时μ0和μ1差距又大乘积才会达到最大。所以最大化类间方差本质上就是在找一个“两类都有足够多的像素且两类均值分得最开”的平衡点。这和最小化类内方差是等价的因为总方差等于类内方差加类间方差总方差是固定的类间方差最大自然类内方差最小。实际写代码的时候用类间方差公式就行不用去算类内方差省事。2.3 Halcon里灰度直方图的数据结构在Halcon里获取灰度直方图最常用的算子是gray_histo它的签名是gray_histo(Regions, Image : : : AbsoluteHisto, RelativeHisto)AbsoluteHisto是绝对直方图就是每个灰度级的像素个数对应上面的n_iRelativeHisto是相对直方图就是每个灰度级的概率对应上面的p_i。这两个都是tuple类型长度等于图像的灰度级数8位图就是256个元素。这里有个细节要注意gray_histo的第一个参数是Regions如果你传的是整幅图像可以用gen_region或者直接用Image的domain。但更常见的做法是先gen_rectangle1画一个ROI然后reduce_domain把图像裁到ROI里再对裁过的图像做gray_histo。这样统计出来的直方图只包含ROI内的像素对于有明确检测区域的场景非常有用能避免背景干扰。还有一个算子叫gray_histo_abs它直接返回绝对直方图不需要指定区域统计的是整幅图像。但如果你要限定ROI还是得用gray_histo配合reduce_domain。另外gray_histo_range可以指定灰度范围适合已知灰度分布区间的场景能减少计算量。3. Halcon实现Otsu的完整代码方案3.1 HDevelop版本从直方图到阈值先上HDevelop的完整实现这段代码可以直接复制到HDevelop里跑。我把它封装成了一个procedure输入是图像和ROI区域输出是Otsu阈值。* 输入Image输入图像Region统计区域 * 输出OtsuThresholdOtsu计算出的阈值 procedure OtsuThreshold(Image, Region : : : OtsuThreshold) * 第一步把图像裁到ROI区域 reduce_domain(Image, Region, ImageReduced) * 第二步计算相对直方图 gray_histo(Region, ImageReduced, AbsoluteHisto, RelativeHisto) * 第三步获取灰度级数 tuple_length(RelativeHisto, NumLevels) * 第四步初始化变量 * 累积概率ω0累积均值μ0的分子部分 * 整幅图均值μT Omega0 : 0.0 Mu0Sum : 0.0 MuT : 0.0 * 先算整幅图均值μT for Index : 0 to NumLevels - 1 by 1 MuT : MuT Index * RelativeHisto[Index] endfor * 第五步遍历所有阈值找最大类间方差 MaxSigmaB : -1.0 BestThreshold : 0 for T : 0 to NumLevels - 2 by 1 * 累积ω0和μ0的分子 Omega0 : Omega0 RelativeHisto[T] Mu0Sum : Mu0Sum T * RelativeHisto[T] * 如果ω0为0或为1跳过避免除零 if (Omega0 1e-10 or Omega0 1.0 - 1e-10) continue endif * 计算μ0和μ1 Mu0 : Mu0Sum / Omega0 Omega1 : 1.0 - Omega0 Mu1 : (MuT - Mu0Sum) / Omega1 * 计算类间方差 SigmaB : Omega0 * Omega1 * (Mu0 - Mu1) * (Mu0 - Mu1) * 更新最大值 if (SigmaB MaxSigmaB) MaxSigmaB : SigmaB BestThreshold : T endif endfor OtsuThreshold : BestThreshold return() endprocedure这段代码的逻辑很直白先算整幅图的均值μT然后从灰度0开始逐个累积ω0和μ0的分子每累积一次就算一次类间方差记录最大值对应的阈值。注意循环只到NumLevels - 2因为阈值t必须保证C1至少有一个灰度级否则ω1为0公式没意义。3.2 关键参数与边界处理上面代码里有几个地方需要特别说明。第一个是1e-10这个容差用来判断ω0是否接近0或1。理论上ω0不会精确等于0或1因为直方图里总有像素但浮点运算可能有微小误差加个容差更稳。第二个是MaxSigmaB初始化为-1.0因为类间方差恒为非负-1保证第一个有效阈值一定会被记录。还有一个容易忽略的点gray_histo返回的RelativeHisto长度是固定的对于8位图是256对于16位图是65536。如果你的图像是16位的这个循环要跑65536次在HDevelop里会明显变慢。这时候有两个优化方向一是用gray_histo_range限定灰度范围只统计实际有像素的区间二是先把图像scale_image或者convert_image_type降到8位牺牲一点精度换速度。具体选哪个看你的检测精度要求。3.3 C接口版本嵌入实际项目HDevelop跑通之后真正上产线通常要用C或者C#调用Halcon库。下面给一个C版本的实现用的是HalconCpp接口。假设你已经有了HImage对象和HRegion对象。#include HalconCpp.h using namespace HalconCpp; int ComputeOtsuThreshold(const HImage image, const HRegion region) { // 裁到ROI HImage imageReduced image.ReduceDomain(region); // 计算直方图 HTuple absHisto, relHisto; imageReduced.GrayHisto(region, absHisto, relHisto); // 获取灰度级数 int numLevels relHisto.Length(); // 计算整幅图均值 double muT 0.0; for (int i 0; i numLevels; i) { muT i * relHisto[i].D(); } // 遍历阈值 double omega0 0.0; double mu0Sum 0.0; double maxSigmaB -1.0; int bestThreshold 0; for (int t 0; t numLevels - 1; t) { omega0 relHisto[t].D(); mu0Sum t * relHisto[t].D(); if (omega0 1e-10 || omega0 1.0 - 1e-10) continue; double mu0 mu0Sum / omega0; double omega1 1.0 - omega0; double mu1 (muT - mu0Sum) / omega1; double sigmaB omega0 * omega1 * (mu0 - mu1) * (mu0 - mu1); if (sigmaB maxSigmaB) { maxSigmaB sigmaB; bestThreshold t; } } return bestThreshold; }C版本和HDevelop版本逻辑完全一致只是语法不同。这里用relHisto[t].D()把HTuple元素转成double注意HTuple的索引是从0开始的和HDevelop一致。实际项目里我会把这个函数放在一个工具类里加上异常处理和日志输出方便排查问题。3.4 算完阈值之后怎么用拿到Otsu阈值之后下一步就是做二值化。Halcon里用threshold算子threshold(ImageReduced, RegionBinary, OtsuThreshold, 255)这里假设是8位图前景是暗的背景是亮的所以阈值下限是OtsuThreshold上限是255。如果前景是亮的就反过来写0, OtsuThreshold。但实际场景往往没这么简单Otsu算出来的阈值可能偏大或偏小这时候可以加一个偏移量threshold(ImageReduced, RegionBinary, OtsuThreshold Offset, 255)Offset根据实际效果调一般正负10到20之间。我一般会先用Otsu算一个基准值然后在界面上留一个滑块让操作员微调这样既保留了自动化的便利又给了人工干预的空间。4. 实操中的性能优化与避坑经验4.1 直方图统计的ROI策略前面反复提到ROI这里展开说一下为什么它这么重要。Otsu的统计是基于直方图的如果整幅图里背景占了90%前景只占10%那直方图里背景的峰值会非常高Otsu算出来的阈值会偏向背景导致前景被切掉一部分。解决办法就是限定ROI只统计你真正关心的区域。ROI怎么画如果是规则工件用gen_rectangle1画个矩形就行如果是圆形工件用gen_circle如果是任意形状可以用gen_region_polygon或者从CAD导入。画好之后用reduce_domain裁图像再传给gray_histo。注意gray_histo的第一个参数也要传ROI区域不能传整幅图否则统计范围就不对了。还有一个进阶技巧如果ROI内还有干扰比如工件上有反光点可以在统计直方图之前先做一次中值滤波或者高斯滤波把噪点压掉。Halcon里用median_image或者gauss_filter窗口大小根据噪点尺寸选一般3x3到7x7。滤波会稍微模糊边缘但对直方图统计影响不大反而能让Otsu阈值更稳。4.2 多通道图像的处理Halcon里彩色图像通常是三通道的gray_histo不能直接对三通道图像用会报错。这时候有两个选择一是用rgb1_to_gray或者rgb3_to_gray转成灰度图再算二是对每个通道分别算Otsu阈值然后取平均或者取某个通道的结果。我一般推荐第一种因为转灰度之后信息虽然损失了一些但Otsu本来就是基于灰度分布的转灰度是标准做法。rgb1_to_gray的权重是0.299R0.587G0.114B符合人眼感知rgb3_to_gray用的是另一种权重。如果工件颜色对比明显比如红色工件在灰色背景上那可以只取R通道或者用access_channel取某个通道效果可能比转灰度更好。4.3 计算速度的实测对比我在一台工控机上做过实测图像是2448x2048的8位灰度图ROI大概占整幅图的四分之一。用HDevelop的binary_threshold算一次大概0.8毫秒用上面手写的Otsu procedure算一次大概2.3毫秒。差距主要在两个地方一是binary_threshold是Halcon底层优化的用了查表和并行二是手写版本在HDevelop里循环256次每次都有tuple索引和浮点运算开销。如果对速度要求高有几个优化方向。第一把循环里的RelativeHisto[T]提前取出来存到一个局部数组里避免每次索引tuple。第二如果灰度分布集中在某个区间用gray_histo_range限定范围循环次数能从256降到几十。第三把HDevelop代码导出成C编译优化之后速度能提升不少实测能到1.2毫秒左右基本够用。第四如果实在要求极致速度可以用Halcon的binary_threshold做粗定位只在它结果附近正负10个灰度级范围内用手写Otsu精算兼顾速度和可控性。4.4 常见问题速查表问题现象可能原因排查方法解决办法阈值算出来是0或255图像几乎全黑或全白直方图集中在一个灰度用min_max_gray看灰度范围检查光源和曝光或先做灰度拉伸阈值偏向背景前景面积太小直方图被背景主导打印直方图看峰值位置缩小ROI或对前景区域加权不同图像阈值波动大光照不稳定直方图整体漂移连续采几帧看直方图变化加光源稳定器或先做emphasize增强对比计算耗时超过5毫秒图像太大或灰度级太多用count_seconds计时缩小ROI降位深或改用C版本二值化后噪点多阈值附近灰度像素多噪声被放大看二值化结果先滤波或阈值加偏移或做形态学开运算16位图循环太慢65536次循环看图像位深用gray_histo_range限定范围或转8位这个表是我自己调试时整理的基本覆盖了八成以上的问题。遇到新问题的时候先看直方图再看ROI最后看位深和光照按这个顺序排查效率最高。4.5 几个容易踩的坑第一个坑是gray_histo的Regions参数和Image参数不一致。比如你传了整幅图作为Image但Regions只传了ROI这时候Halcon会报错或者结果不对。正确做法是先用reduce_domain把图像裁到ROI然后gray_histo的两个参数都传ROI。第二个坑是浮点图像的直方图。如果图像是real类型gray_histo会报错必须先转成byte或者uint2。转的时候注意scale_image的范围别把有用信息截断了。第三个坑是阈值边界。循环里如果写成for T : 0 to NumLevels - 1当T取到最后一个灰度级时ω1为0除零会报错。所以一定要写到NumLevels - 2并且在循环里加ω0和ω1的容差判断。第四个坑是HDevelop和C的tuple索引差异。HDevelop里tuple索引从0开始C里HTuple索引也从0开始但如果你用HTuple::operator[]返回的是HTupleElement要调.D()转double。这个细节在混合编程的时候容易搞混。5. 从Otsu延伸出去的实用技巧5.1 自适应Otsu分块计算如果图像光照不均匀整幅图用一个Otsu阈值效果不好这时候可以把图像分成若干小块每块单独算Otsu阈值然后对每个像素用其所在块的阈值做二值化。这就是自适应Otsu的思路类似local_threshold但更可控。具体做法把图像用crop_domain或者tile_images分成MxN块每块调一次Otsu procedure把阈值存到一个矩阵里。然后遍历每个像素根据坐标找到对应的块阈值做判断。Halcon里可以用get_grayval逐像素读但这样太慢更好的做法是用threshold对每块分别二值化然后用concat_obj或者union2把结果拼起来。分块大小一般选64x64到128x128太小了每块像素不够直方图统计不准太大了又失去了自适应的意义。5.2 多级Otsu不止分两类标准Otsu是分两类但有些场景需要分三类甚至更多比如工件、背景、缺陷三种区域。多级Otsu就是把类间方差的概念扩展到多个阈值找一组阈值使多类的类间方差最大。实现上可以用动态规划或者穷举但计算量比两级大很多。实际项目里我一般用两级Otsu先分前景背景然后在前景里再跑一次Otsu分缺陷和正常区域这样计算量可控效果也够用。5.3 结合形态学做后处理Otsu算出来的二值化结果往往有毛刺和孤立噪点这时候用形态学处理一下。opening_circle或者opening_rectangle1去掉小噪点closing_circle填补小孔洞。结构元素大小根据缺陷尺寸选一般比最小缺陷小一半左右。注意形态学会改变区域边界如果后续要测尺寸得在形态学之前先做测量或者用shape_trans把区域转成标准形状再测。5.4 阈值稳定性监控产线上最怕的是阈值突然跳变导致误判。我一般会在程序里加一个监控每次算完Otsu阈值和上一次的阈值比较如果变化超过某个范围比如20个灰度级就报警或者触发重新标定。同时把每次的阈值和直方图峰值位置记到日志里方便追溯。Halcon里可以用write_tuple把阈值写到文件或者用open_file写日志。这个习惯帮我提前发现过好几次光源老化的问题。5.5 和Halcon自带算子的对比选择最后说一下什么时候用自写Otsu什么时候用binary_threshold。如果项目用的是较新版本的Halcon图像是标准8位灰度光照稳定ROI明确那直接用binary_threshold的max_separability模式最省事速度和稳定性都是最优的。但如果遇到下面几种情况就值得自己写老版本Halcon算子行为异常需要限定ROI但binary_threshold不支持需要记录中间变量做调试需要做多级或者自适应扩展需要把算法移植到不支持Halcon的平台上。我自己在产线项目里的策略是先用binary_threshold快速验证如果效果和速度都满足就保留如果发现阈值不稳定或者需要定制再换成手写版本把参数暴露出来给操作员调。我个人在实际操作中的体会是Otsu这个算法看起来简单但真正写好、用稳关键不在算法本身而在对图像的理解和对场景的把控。直方图长什么样、ROI怎么画、光照怎么控、阈值怎么监控这些才是决定成败的地方。代码只是最后一步前面的准备工作做足了代码自然就稳了。