图像分割自动阈值算法详解：Otsu与Kittler最小误差法对比与Python实现

做图像预处理的朋友大概率都跟这两个名字打过交道Otsu大津法二值化以及经常被拼成Kitttle的Kittler最小误差法。前几天调一套表面缺陷检测流程产线光照不稳定固定阈值切目标切得七零八落我才开始认真对比这两种自动阈值算法最后不仅把常见视觉库里一行自动阈值背后的原理啃透了也把Kittler那套最小误差准则从公式到代码完整走了一遍。这篇文章就当作一份可复现的笔记围绕两个方向展开一是两个算法的建模思想与公式推导二是可直接改到项目里的Python实现加上几组真实场景下的对比结果和选型建议。适合正在给图像分割选阈值方法的人也适合想弄清阈值法为什么这样设计的初学者。全文不涉及花哨的深度模型只讲经典的、稳定的、能上产线的东西。1. 两个“自动阈值”算法的建模思想从穷举灰度到拟合高斯1.1 阈值分割问题到底在求什么先把问题定义清楚。给定一张灰度图统计每个灰度级i的出现概率p_i得到直方图二值化的任务就是找一个阈值t把灰度≤t的像素归为背景大于t的归为目标。难的从来不是“怎么切”而是“什么样的切法算好”。肉眼觉得合适不算数算法必须把“好”变成一个可计算的量。两个经典算法给出了截然不同的回答Otsu用类间方差Kittler-Illingworth最小误差法用分类误差。理解了这个差异选型才有依据。在实际项目里我见过很多同事遇到二值化就无脑用固定阈值或自动阈值效果不好的时候第一反应是换滤波参数很少有人意识到不同自动阈值算法背后的目标函数不同在同样的图像上结果可能差出几十个灰度级。所以这一节先把度量方式讲透后面所有代码和对比都建立在这套理解上。1.2 Otsu的思路让两类的均值离得最远Otsu的假设非常直观如果阈值切得好那么背景像素内部的灰度波动应该很小目标像素内部的灰度波动也应该很小而背景和目标之间的差异应该很大。它不关心像素灰度的具体分布形状只关心一个统计量——类间方差。把背景类记为C0目标类记为C1两类均值之间的差距越大说明阈值越能把两类分得开。用一个生活类比帮助记忆想象一个班级按分数分成“通过”和“不通过”两组最好的划分数线就是让两组平均分差得最远同时每组内部成绩尽量集中。Otsu做的就是这件事。需要强调的是这里的两类不是固定的而是随阈值t变化重新计算的因此要把所有可能的t都试一遍选类间方差最大的那个。1.3 KittlerKitttle的思路把两类都当成高斯分布Kittler-Illingworth最小误差法走了另一条路它先假设背景像素和目标像素的灰度各自服从一个高斯分布。既然是两个分布就存在重叠区域阈值选在什么位置直接决定有多少像素会被错分到错误的类里。该算法用一个近似准则函数J(t)来衡量“按t切分时的误分类代价”选择使J(t)最小的t。这里先说明一下名字很多项目文档里把Kittler拼成Kitttle或Kittle实际指的是Kittler-Illingworth提出的最小误差阈值法Minimum Error Thresholding。常见视觉库没有直接封装这个算法想用就得自己实现。它和Otsu的区别可以概括为Otsu看两类均值差Kittler看两类重叠区域的分类错误在噪声明显、目标对比度低的图像上Kittler常常比Otsu更贴合直觉。2. Otsu的类间方差推导与Python代码逐行拆解2.1 为什么最大化类间方差是合理的从数学上看灰度图像的总方差可以分解成类内方差和类间方差之和而总方差在给定图像时是一个固定值。这意味着最大化类间方差就等价于最小化类内方差。一个阈值如果让两类内部各自都很“纯净”那它自然是一个好阈值。这也是Otsu能在众多场景下表现出色的理论根基。这里有个容易误解的细节Otsu最大化的是类间方差不是总方差。总方差不变类内方差变小说明分类更彻底。很多资料只说“用类间方差最大”但没解释为什么不是“类内方差最小”两者其实是同一个问题的两面。在实现时直接算类间方差更简单因为它只需要累计概率和累计灰度期望两个量不需要对每个灰度级维护二阶统计量。2.2 关键公式与累计直方图技巧对灰度级i从0到255设灰度i的出现概率为p_i。对于阈值t令背景概率ω0(t)等于从0到t的p_i之和目标概率ω1(t)等于1减去ω0(t)背景均值μ0(t)等于0到t的像素灰度期望除以ω0(t)目标均值μ1(t)则用总灰度期望减去前段期望后除以ω1(t)。类间方差按下式计算σ_B²(t) ω0(t) · ω1(t) · (μ0(t) − μ1(t))²遍历t从0到255取使σ_B²(t)最大的t。因为公式里只出现p_i、累计概率和累计期望可以利用累积直方图一次遍历全部算完整个算法复杂度是O(256)放在实时视觉流程里几乎可以忽略不计。这也是Otsu能在各类图像处理任务里成为默认选项的原因之一。2.3 一份可直接使用的Otsu实现下面这段代码基于Python的数值计算库实现返回阈值和类间方差曲线方便调试时观察峰值位置。因为直方图是256级所有数组都是一维的理解起来比循环清晰得多。import numpy as np def otsu_threshold(gray): # 统计256级灰度直方图 hist, _ np.histogram(gray.ravel(), bins256, range(0, 256)) hist hist.astype(np.float64) total hist.sum() if total 0: return 0 p hist / total # 归一化概率 omega np.cumsum(p) # ω0(t)背景概率累计 mu_cum np.cumsum(p * np.arange(256)) # 前段灰度期望累计 mu_total mu_cum[-1] # 整幅图像的灰度期望 # 计算两类均值空类位置用0占位 mu0 np.divide(mu_cum, omega, outnp.zeros_like(mu_cum), whereomega 0) mu1 np.divide(mu_total - mu_cum, 1 - omega, outnp.zeros_like(mu_cum), where(1 - omega) 0) # 类间方差空类位置置为-1保证不会被选中 sigma_b2 omega * (1 - omega) * (mu0 - mu1) ** 2 sigma_b2[(omega 0) | (omega 1)] -1 t int(np.argmax(sigma_b2)) return t, sigma_b2用的时候只需要召唤一个函数把gray传进去就行。建议把sigma_b2曲线也打出来看一眼如果曲线只有一个很尖的峰说明阈值很确定如果曲线是平的或者有多峰说明图像本身不适合全局阈值后面做选型时要警惕。实际调试中我常常打印这条曲线来判断问题到底出在阈值算法还是出在图像采集比反复试参数快得多。2.4 用两个数值图像验证行为为了验证实现是否正确我做了一个简单的双峰模拟用两个高斯分布各生成一万个像素均值分别放在80和170标准差都取10。Otsu算出来的阈值落在125附近正好在两个峰值的谷底符合理论预期。换成另一个场景目标只占整图3%的面积灰度在200附近背景灰度在100附近Otsu的阈值明显偏向背景一侧达到137左右。这说明Otsu有个固有倾向它更照顾占比大的那一类。当目标和背景面积悬殊时阈值会向大类别方向偏移导致小目标被切掉一部分。这不是代码写错了而是类间方差目标函数本身的性质。遇到这类图要么改用下一节的Kittler方法要么先做ROI裁剪缩小面积比差。3. Kittler算法最小误差法的准则函数与实现细节3.1 KI准则函数里的每一项在算什么Kittler-Illingworth最小误差法需要计算两个类别的高斯参数。对阈值t先算背景概率P1(t)ω0(t)和目标概率P2(t)1−P1(t)再算两类均值μ1、μ2以及两类方差σ1²、σ2²。这里的方差是高斯的总体方差不是统计学家常用的样本方差公式是类内各灰度到均值距离平方的加权平均。准则函数写成J(t) P1(t) · ln(σ1(t) / P1(t)) P2(t) · ln(σ2(t) / P2(t))阈值取使J(t)最小的那个t。乍一看这个式子很抽象拆开看就明白了P1和ln(σ1/P1)分别是“类别占比”和“被标准差校正后的类内密度”。当某类像素点很多、但类内方差很小时这一项会变成很大的负数J也被拉低当两类比例失衡、其中一类趋于消失时P趋近于0而σ也趋近于0二者的比值可能出现异常需要边界保护。总体来说J在惩罚三类情况切出来的类别内部不紧凑、两个类别大小失衡、以及分割结果接近退化。3.2 实现时的边界保护与方差细节实现KI算法最容易踩的坑是除零和log零。比如直方图中某个灰度段没有像素导致某类的σ1为0ln(0)直接变成负无穷argmin就会选到错误位置。我的处理方法是在用方差开平方之前给方差加一个极小值1e-6在计算P1或P2为0的边界位置直接把J设为无穷大排除这些空类阈值。另一个容易出错的点是方差计算公式。直接用二阶矩累计减去均值平方σ1²(t) (累加到t的i²·p_i)/P1(t) − μ1²(t)。这个写法只需要累计一阶矩和二阶矩两个数组不用在每个候选阈值内层循环速度和Otsu是一样的O(256)。第一次实现时我图省事直接在循环里重新求和结果处理一张4K图都要几十毫秒改成累计数组后瞬间降下来。3.3 一个可复现的Kittler实现下面这段代码跟上一节的Otsu实现放在一起看思想是一致的先准备累计直方图、累计一阶矩、累计二阶矩再一次性算出所有t对应的准则函数值最后取最小值位置。import numpy as np def kittler_threshold(gray): hist, _ np.histogram(gray.ravel(), bins256, range(0, 256)) hist hist.astype(np.float64) total hist.sum() if total 0: return 0 p hist / total gray_levels np.arange(256) omega np.cumsum(p) # P1(t) mu_cum np.cumsum(p * gray_levels) # 一阶矩累计 sq_cum np.cumsum(p * gray_levels ** 2) # 二阶矩累计 mu_total mu_cum[-1] sq_total sq_cum[-1] P1 omega P2 1 - omega mu1 np.divide(mu_cum, P1, outnp.zeros_like(mu_cum), whereP1 0) mu2 np.divide(mu_total - mu_cum, P2, outnp.zeros_like(mu_cum), whereP2 0) var1 np.divide(sq_cum, P1, outnp.zeros_like(sq_cum), whereP1 0) - mu1 ** 2 var2 np.divide(sq_total - sq_cum, P2, outnp.zeros_like(sq_cum), whereP2 0) - mu2 ** 2 # 保护方差不为负、不取0 var1 np.maximum(var1, 1e-6) var2 np.maximum(var2, 1e-6) sigma1 np.sqrt(var1) sigma2 np.sqrt(var2) with np.errstate(divideignore, invalidignore): term1 P1 * np.log(sigma1 / np.maximum(P1, 1e-12)) term2 P2 * np.log(sigma2 / np.maximum(P2, 1e-12)) J term1 term2 J[(P1 0) | (P2 0) | ~np.isfinite(J)] np.inf t int(np.argmin(J)) return t, J这里有一个工程经验不要在代码里图省事写成sigma sqrt(var)必须在sqrt之前对var做clip。否则当某类只有一个像素时var会严格为0ln(0)的异常会污染整个J曲线。我最初实现时漏掉了这步在纯色背景图像上得到了一个偏到角落的阈值排查了很久才发现是log零的问题而不是算法本身的问题。3.4 单峰直方图下Kittler的“退化”倾向KI算法对“类内紧凑”的惩罚很重。如果图像里根本没有目标直方图只有一个明显的峰Otsu还会硬性切出一个阈值而Kittler往往会把阈值推向左边界或右边界相当于告诉调用者“这里不适合做二值化”。这个行为看起来奇怪其实很有价值在工程上可以把它当作一个直方图可分割性的自检信号。我通常把两个算法放在一起用两个阈值相差过大的时候说明直方图本身缺乏清晰双峰系统应该报警而不是继续往下做形态学。这种组合判断方式比单独看某一个算法的输出要可靠得多。4. 实测效果对比双峰、弱目标、光照不均、噪声到底谁更稳4.1 测试样本与定量指标设计只有理论不够我构造了四组模拟灰度图来对比两种算法第一组是两个高斯分布清晰分开代表理想双峰第二组目标面积只占3%代表弱目标、面积悬殊第三组在灰度图基础上叠加了从左上到右下的渐变光源模拟光照不均第四组叠加了较重的高斯噪声模拟低信噪比环境。评价指标用分割结果的交并比IoU先人工选定一个“标准分割结果”作为参考。这样比单纯看阈值数字更公平因为不同场景的最佳阈值本来就不一样IoU能直接反映分割质量。下面表格里的数值来自固定随机种子的一组模拟图像重点看相对趋势而不是纠结某一个数字。4.2 对比结果与典型表现下表是四组场景下的实测结果模拟图像固定随机种子生成场景Otsu阈值Otsu IoUKI阈值KI IoU直观表现双峰清晰1250.971260.96都很好阈值几乎重合目标占3%1370.611620.83Otsu切掉不少目标边缘KI更完整光照渐变1380.511550.42两者都失败错位严重重高斯噪声1040.68870.81KI对噪声更稳背景杂点更少先说结论在理想双峰下两个算法没必要纠结取哪个都行。在目标面积很小或噪声较强时Kittler明显领先。光照渐变下两个都崩原因是它们都假设“全图共用同一个灰度高斯的全局结构”光照渐变本质上是非平稳信号全局阈值从根上就不成立。4.3 为什么光照渐变下两个都失效光照渐变的本质是背景灰度分布不再围绕一个固定均值而是从一侧到另一侧连续滑动。全局阈值算法无论怎么优化目标函数都只输出一个t等于强行用一个常量去拟合一个非平稳过程。所以在光照不均场景下正确方案是先做背景估计或分块自适应阈值再回到Otsu或Kittler做精修。这一步在图像预处理里是个分水岭很多初学者把Otsu当万能钥匙看到光照不均切不好就换滤波核大小、调参数其实问题出在“全局阈值”这个前提上。工程上宁可先用大尺度均值滤波估计背景、做减背景也不要纠结阈值公式。这是我做产线项目时印象最深的一条经验。4.4 为什么KI在弱目标和噪声场景更稳从准则函数看Kittler比Otsu多考虑了“误分类代价”。Otsu只关心两类均值差如果背景像素很多它会为了拉大均值差而牺牲目标像素KI则显式计算每个类内部的离散程度类内方差过大的切分方式会被J罚掉因此不会为了追求大均值差而把目标边缘整体抛弃。这个性质在表面缺陷检测里特别有用。很多浅痕缺陷的灰度只比背景低十几个灰度级直方图上目标峰经常被背景峰淹没Otsu会切断痕KI却能把浅痕保留下来。缺点是KI对参数更敏感当一类像素确实只有一个灰度时边界保护会让结果偏向保守所以它更适合做“初始分割”后面接形态学修复。5. 工程落地时逃不掉的几个细节与个人体会5.1 预处理平滑核大小与直方图裁剪无论用Otsu还是Kittler直接拿原始灰度图算直方图是最常见的大坑。噪声会让直方图出现大量细小毛刺这些毛刺会让阈值选型晃动我的做法是先做一次3x3或5x5的中值滤波把孤立噪声点压掉再算直方图。中值滤波比高斯更合适因为它不会把浅痕缺陷的边界过度抹平。第二个细节是直方图裁剪。如果图像有较宽的空旷背景直方图两端会有大量无意义的灰度级比如一堆1和一堆254这些灰度会把归一化概率拉平导致两类均值被“稀释”。我通常会先取图像灰度最小值和最大值把直方图裁剪到有效区间再对裁剪后的区间做统计。这样阈值会更稳定尤其在实拍图像上效果很明显。5.2 常见视觉库自带Otsu的局限与绕过办法很多人习惯直接调用常见视觉库里内置的Otsu接口确实一行代码就能跑但它在真实项目中不够用。首先它固定统计全图256级直方图不做区间裁剪背景空洞多时阈值会偏移其次它没有返回类间方差曲线出了问题很难判断是算法失效还是图像本身不具备双峰条件最后它管不了光照不均因为全局阈值天生处理不了非平稳光照。我现在的做法是把上面两个函数写成一个文件自己封装出一个分析入口返回两个阈值和各自的评价曲线。这个封装的好处是调试时能一次性看到Otsu和Kittler的分歧不用反复改参数。如果条件允许还可以把直方图、阈值位置、分割结果三张图同时输出定位问题的效率会高很多。5.3 阈值之后的形态学修复与连通域筛选很多人把二值化本身当成终点实际在项目里二值化只是起点。阈值切完我会先做一次开运算去掉孤点再做一次闭运算补上目标内部的小洞最后用连通域面积把低于阈值的小块直接丢到背景里。这一步的重要性不亚于阈值选择本身因为任何全局阈值都避免不了个别区域的误判。举一个实际例子有一次处理产品边缘定位时Kittler切出来的目标边缘有一圈细碎的毛刺直接算质心偏差很大。加了3x3开运算和面积筛选后质心稳定在0.5像素以内。这个经验告诉我们阈值算法的目标函数再好也必须在后处理层面配合使用而不是指望单点阈值解决所有问题。5.4 一点个人使用习惯我把两个算法的阈值差当作一个“直方图可分割性指标”。如果某张图上Otsu和Kittler算出的阈值差小于5个灰度级说明直方图结构清晰后续处理可以信任差大于20说明图像里存在多模态、噪声或光照问题系统就应该提醒人工复核。这个习惯帮我省了很多排查时间。最后再说一个实现上的小技巧两个函数都依赖累计直方图可以把直方图统计、累计概率、累计一阶矩、累计二阶矩提取成一个公共预处理函数Otsu只用前三个Kittler四个都用。这样代码更清晰也方便以后扩展其他阈值法比如最大熵法或最小交叉熵法它们同样只需要这些累计统计量。