做图像分割时最朴素也最常用的自动阈值方法就是大津法也就是 OTSU 算法。我在处理文档扫描件、工业零件检测图像时第一步基本都会用 OTSU 把前景和背景分开Python 里配合 OpenCV-Python 一行就能跑出阈值确实方便。但如果你直接把一维 OTSU 用在噪声偏大、背景不匀的图上很容易出现目标被吞掉、分割区域粘连的情况。后来我把算法升级成二维 OTSU在计算像素灰度阈值的同时把邻域平均灰度也纳入进来等于同时看“单个点”和“周围一片”的信息抗噪能力明显提升。这篇文章就结合我的 Python 手写实现和 OpenCV-Python 环境下的真实调用过程把二维 OTSU 的原理、公式、代码、优化和避坑经验完整讲一遍。2. 项目背景为什么我盯上二维 OTSU2.1 一维大津法的基础与局限先快速回顾一下一维大津法。OTSU 由日本学者大津展之在 1979 年提出核心思想是按灰度阈值把图像分成背景和目标两类然后计算两类之间的类间方差让这个方差最大的灰度值就是最佳阈值。类间方差越大说明两类区分越明显分割效果通常也越好。用 Python 实现的时候一维 OTSU 的计算过程其实很直观统计灰度直方图然后从 0 到 255 遍历阈值 t对每个 t 计算背景概率、目标概率、两类各自的灰度均值最后套用类间方差公式。OpenCV-Python 里更省事直接调 cv2.threshold 的 THRESH_OTSU 标志位就能拿到阈值不需要自己写。但一维 OTSU 有个绕不开的短板它只用单个像素的灰度值做判断完全没有考虑像素之间的空间关系。只要图像里混入噪声或者在光照不均匀的背景下单个像素灰度很容易“误判”。我实测过一张带高斯噪声的文档图一维 OTSU 分完以后文字区域里面出现大量白色空洞背景里反而残留不少黑色噪点。这种情况下面准确率再高也只是在单点特征上做文章救不回来。2.2 二维 OTSU 能解决什么问题二维 OTSU 的出发点很简单把每个像素的信息从“当前灰度值”扩展成“当前灰度值 邻域平均灰度”的二元组。为什么加邻域平均因为真正的目标和背景在空间上通常是连片的目标内部和背景内部的灰度变化相对平缓噪声点往往是孤立的。把邻域信息带进来以后孤立噪声在二维特征上就会偏离主体区域更容易被识别出来、被滤除。我自己的体会是二维 OTSU 在处理三类图时特别有价值一是低对比度文档图像二是工业表面缺陷检测中的纹理噪声干扰三是医学影像里目标和背景灰度范围重叠较多的场景。它比一维版本多了一个维度的上下文信息误分割率明显下降代价是计算量成倍增加原理理解起来也稍微绕一些。3. 二维 OTSU 的原理与关键推导3.1 把图像信息拓展成二维特征一维 OTSU 使用直方图统计的是“灰度值 i 出现了多少次”二维 OTSU 统计的是“灰度值 i 且邻域平均灰度为 j 的像素有多少个”。这个统计结果可以画成一个 256 x 256 的二维直方图。邻域平均灰度怎么算常见做法是用一个固定大小的窗口对图像做均值滤波。窗口尺寸通常是 3x3 或 5x5我常用 5x5抗噪能力更强后面细说。滤波之后的图像每个像素值代表原图像素周围一片的平均灰度可以理解为一种“局部背景估计”。于是每个像素都拿到了一个特征向量 (i, j)i 是原始灰度j 是局部均值。理想情况下如果把每个像素画在 i-j 坐标系上你会发现大量像素点集中在对角线附近因为正常情况下像素值和它周围的平均值差异不大。而噪声点、边缘点由于局部突变会明显偏离对角线。二维 OTSU 要做的就是在这个二维坐标系里画一条“L 形”分界线把背景和目标分开。3.2 二维直方图与区域划分假设灰度级数仍为 L256二维直方图用矩阵 H[i][j] 表示元素值为满足“灰度i 且邻域均值j”的像素个数。把 H 除以总像素数得到归一化概率 p(i,j)。现在设置一个阈值向量 (s, t)其中 s 是灰度阈值t 是邻域均值阈值。以 (s, t) 为分界点二维直方图被对角线方向的十字结构划分成四个区域区域 Ai ≤ s 且 j ≤ t通常对应背景。区域 Bi s 且 j t通常对应目标。区域 Ci ≤ s 且 j t多是边缘或噪声。区域 Di s 且 j ≤ t同样是边缘或噪声。大多数二维 OTSU 实现把 C 和 D 当成干扰区域忽略掉因为真正目标和背景的主体都紧贴对角线分布。计算的时候只统计 A 和 B 两个区域的概率分布。这样做有个隐含前提噪声和边缘像素数量占比不大忽略后对全局中心的影响可以接受。3.3 类间方差度量与阈值求解有了区域划分接下来要定义“两类分得够不够开”的度量。一维 OTSU 用灰度均值之差构造类间方差二维 OTSU 则用两个均值向量。记背景区域 A 的概率为 w0目标区域 B 的概率为 w1。背景的灰度均值向量为 u0 (u0i, u0j)目标的灰度均值向量为 u1 (u1i, u1j)。再计算全局均值向量 uT (ui, uj)其中 ui 是全图灰度均值uj 是全图邻域均值。二维 OTSU 的目标函数通常取类间离散度矩阵的迹也就是tr(Sb) w0 * [(u0i - ui)^2 (u0j - uj)^2] w1 * [(u1i - ui)^2 (u1j - uj)^2]遍历所有可能的 (s, t)使 tr(Sb) 最大的那一组 (s, t) 就是最佳阈值向量。这个式子的直观含义是背景中心、目标中心与全局中心的加权距离越大说明两类分离度越高分割效果越好。这里有一个细节需要注意由于 C、D 区域被忽略w0 w1 不一定等于 1。计算目标区域 B 的均值时要单独对 B 区域内的像素求和再除以 w1不能拿“全局累计减去 A 区累计”直接当 B 区累计用否则会把 C、D 区域的贡献混进来导致阈值偏差。这个坑我在早期实现里踩过后面代码部分会体现正确的处理方式。4. Python 手写实现全流程4.1 环境准备与图像预处理开发环境我建议直接用 Python 3.8 以上版本配上 numpy 和 opencv-python。如果还没装 OpenCV-Python在终端执行 pip install opencv-python 即可。numpy 一般装 OpenCV 时也会自动带进来但为了保险可以显式装一下。准备一张灰度测试图尺寸别太大先用 300 x 300 左右比较好跑。如果是彩色图先转成灰度import cv2 import numpy as np img cv2.imread(test.jpg) gray cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)这里要特别注意二维 OTSU 的输入图像必须是 8 位单通道图也就是灰度值范围 0 到 255。如果你拿 16 位图直接算直方图维度就不匹配了。4.2 邻域均值和二维直方图构建邻域均值我用 OpenCV 的滤波函数计算。它比单纯用卷积核要稳尤其是边界处理可以直接指定 borderType。常用的边界策略是 BORDER_REPLICATE也就是复制边缘像素避免边界处因为补零导致均值偏低。win 5 kernel np.ones((win, win), np.float32) / (win * win) mean cv2.filter2D(gray, -1, kernel, borderTypecv2.BORDER_REPLICATE) mean np.round(mean).astype(np.int32)窗口大小对结果影响很大。3x3 保留细节但抗噪一般5x5 是折中7x7 更平滑但容易把细线条目标模糊掉。我处理文档文字时常用 3x3处理工业零件表面时用 5x5具体要根据图像分辨率来调。接下来构建二维直方图。第一版我写过双循环遍历每个像素600 x 600 的图跑了差不多 1 秒多勉强能接受。但像素更多就太慢了后面会讲优化方案。基础版本代码如下h, w gray.shape hist np.zeros((256, 256), dtypenp.float64) for r in range(h): for c in range(w): i gray[r, c] j mean[r, c] hist[i, j] 1.0 hist / (h * w)这里有个细节mean 数组要转成整数再去索引直方图不能直接拿浮点数当索引。我最初在这里直接用了 mean 的浮点值结果程序报错后来才发现要 round。4.3 积分图加速与完整核心代码如果按最朴素的思路对每个候选 (s, t) 都重新扫描 256 x 256 的直方图区域来求 w0、u0、u1那总复杂度是 O(L^4)。L256 时256 的 4 次方约 43 亿次操作放到任何解释型语言里都不可接受。所以必须引入二维前缀和也就是积分图。积分图的好处是任意矩形区域的像素和都能在 O(1) 时间内取出来。对于二维 OTSU我们需要维护三个积分图概率累积、灰度值乘以概率的累积、邻域均值乘以概率的累积。构建方式可以用嵌套循环也可以直接用 numpy 的 cumsum。我用 numpy 的 cumsum 构建代码简洁还不容易写错def integral_image(mat): return np.cumsum(np.cumsum(mat, axis0), axis1)然后对 hist、i * hist、j * hist 分别求积分图。之后遍历所有 (s, t) 时用积分图快速算出 A 区域和 B 区域的累积量。这里放完整实现def otsu_2d(gray, win5): # 1. 邻域均值 kernel np.ones((win, win), np.float32) / (win * win) mean cv2.filter2D(gray, -1, kernel, borderTypecv2.BORDER_REPLICATE) mean np.round(mean).astype(np.int32) # 2. 二维直方图 hist np.zeros((256, 256), dtypenp.float64) h, w gray.shape for r in range(h): row_gray gray[r] row_mean mean[r] for c in range(w): hist[row_gray[c], row_mean[c]] 1.0 hist / (h * w) # 3. 构建三个积分图 idx np.arange(256, dtypenp.float64) p_integral np.cumsum(np.cumsum(hist, axis0), axis1) ip_integral np.cumsum(np.cumsum(hist * idx.reshape(-1, 1), axis0), axis1) jp_integral np.cumsum(np.cumsum(hist * idx.reshape(1, -1), axis0), axis1) total_p p_integral[255, 255] total_i ip_integral[255, 255] total_j jp_integral[255, 255] ui total_i / total_p uj total_j / total_p # 4. 遍历阈值 max_sigma 0.0 best_s 0 best_t 0 for s in range(256): for t in range(256): # A 区域is, jt p_a p_integral[s, t] if p_a 0: continue i_a ip_integral[s, t] / p_a j_a jp_integral[s, t] / p_a # B 区域is, jt p_b total_p - p_integral[s, 255] - p_integral[255, t] p_integral[s, t] if p_b 0: continue i_b (total_i - ip_integral[s, 255] - ip_integral[255, t] ip_integral[s, t]) / p_b j_b (total_j - jp_integral[s, 255] - jp_integral[255, t] jp_integral[s, t]) / p_b sigma p_a * ((i_a - ui) ** 2 (j_a - uj) ** 2) \ p_b * ((i_b - ui) ** 2 (j_b - uj) ** 2) if sigma max_sigma: max_sigma sigma best_s s best_t t return best_s, best_t注意看代码里 B 区域的取法。用 total 减去“is 的所有 j”、“所有 i 且 jt”再加回左上角区域得到的是严格意义上的 is 且 jt 区域这样就把 C、D 区域排除在外了。我当时第一次实现时直接用 total_p - p_a 当 p_b结果阈值始终偏大后来打印每个区域的概率才发现 B 区域把右下以外的大量噪声点也算进去了。4.4 运行结果示例我拿一张加了高斯噪声的简单二值图测试。原始图是一块黑色背景上有一个白色矩形噪声标准差大约 20。一维 OTSU 返回的阈值大约在 128分割后白色矩形边缘出现很多毛刺背景里也有不少残留噪点。二维 OTSU 用 5x5 窗口跑出来的结果是 s96t128按这两个阈值分割后噪声点大部分被滤掉了矩形边缘干净很多。当然每个图像的结果都不一样这个数据只是参考。你可能会问为什么 s 和 t 不相等因为原图灰度值和邻域均值的分布本来就是两个维度最佳阈值自然不一定落在对角线上。这个特点正是二维 OTSU 比一维灵活的地方。5. 在 OpenCV-Python 中的实际应用与对比5.1 使用 cv2.threshold 的 OTSU 模式OpenCV 自带一维 OTSU调用方式很简单thresh, binary cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU)这里 threshold 的第一个参数传 0 就行因为 THRESH_OTSU 会忽略指定阈值自行搜索最佳阈值。返回的第一个值就是最终阈值第二个值就是分割后的二值图。OpenCV 并没有内置二维 OTSU 的接口所以二维版本必须自己实现然后把算出来的 (s, t) 用于分割。分割时不是简单地把灰度大于 s 的像素设为 255 就完了而是要同时满足灰度阈值和邻域阈值这样才真正用上二维信息。分割策略有几种。最常用的是像素同时满足 gray s 且 mean t才判为目标否则判为背景。用代码写就是s, t otsu_2d(gray, win5) kernel np.ones((5, 5), np.float32) / 25 mean cv2.filter2D(gray, -1, kernel, borderTypecv2.BORDER_REPLICATE) binary np.where((gray s) (mean t), 255, 0).astype(np.uint8)你也可以用另一种方式先分别用两个阈值做二值图再取交集或并集。实际测试下来直接按“双条件”分割效果最稳定。5.2 一维与二维分割效果对比为了直观我做了三组对比。一组是普通文档扫描文字一组是带椒盐噪声的零件图一组是低对比度的医学灰度图。文档文字那张图一维 OTSU 已经表现不错文字边缘偶尔有断裂二维 OTSU 补上了一些断裂处整体差别不算大。但在椒盐噪声的零件图上差距就明显了。一维结果里噪声点大量残留二维结果干净很多。原因就在于椒盐噪声点的灰度值往往极低或极高但邻域均值与主体差别很大二维条件下会被排除掉。低对比度医学灰度图那一组一维 OTSU 很容易把目标区域和背景灰度接近的部分分割错二维 OTSU 因为结合了局部一致性轮廓更完整。不过要注意二维 OTSU 处理低对比度图也不是万能的如果目标和背景的邻域均值分布高度重叠照样会失效。5.3 参数调优与窗口选择窗口 size 是二维 OTSU 最重要的超参数。我做过一个小实验分别用 3x3、5x5、7x7 窗口跑同一张噪声水平不同的图。窗口越小二维直方图越接近对角线结果越像一维 OTSU抗噪能力弱一点但目标细节保留得好。窗口越大邻域均值越平滑抗噪强但目标边缘容易被“抹平”分割出来的边界会往目标内部缩。分辨率较高的图像我一般选 5x5分辨率低或者目标很细的时候选 3x3目标是大块均匀区域时可以上 7x7。还有一个容易被忽略的点邻域均值计算方式。如果直接用 cv2.boxFilter 或 cv2.filter2D 做均值结果一样。但如果用高斯模糊代替均值模糊邻域均值会带有加权意义结果可能略有不同。我测试下来差别不大常规情况就用均值。6. 踩坑记录与工程化建议6.1 直方图构建太慢的问题前面代码里的二维直方图用 Python 双循环构建对 600 x 600 的图像大约要 1.2 秒左右。如果图像到 2000 x 2000这个循环会变成瓶颈耗时十几秒。优化思路有几个。第一是降采样灰度级把 256 级压缩到 64 级甚至 32 级直方图矩阵变小遍历阈值的开销也大幅下降。压缩方式可以用灰度值除以压缩系数取整。第二是空间降采样先缩小图像尺寸计算邻域均值和直方图再把阈值映射回原尺寸。第三是用 Cython、numba 或 C 重写直方图构建部分。我实际工程里最常用的是“灰度级压缩 numba 加速”组合。把灰度压缩到 64 级后阈值遍历次数从 65536 降到 4096速度提升约 16 倍。配合 numba 的 jit 装饰器整个流程从 1 秒多降到几十毫秒应付视频帧里的实时分割不太够但处理单张图片已经很舒服了。6.2 阈值结果不稳定的排查二维 OTSU 偶尔会出现两帧相近图像算出的阈值跳变很大的情况尤其是有大量噪点或图像内容变化剧烈时。排查思路首先看二维直方图是否异常。如果图像中某个灰度区间的像素极少积分图里对应区域概率接近 0可能导致均值计算出现极端值。其次看邻域窗口是否过大窗口大了以后邻域均值过于平滑二值图里小目标容易被吞掉阈值随之漂移。最后看 C、D 区域是不是占比太高。如果噪声和边缘像素占比超过一定比例忽略它们的做法就不再可靠。一个比较实用的稳定化手段是对输入图像先做轻微高斯滤波减少极端噪声再把二维 OTSU 得到的结果与一维 OTSU 结果做加权融合避免在工况差异大的场景里频繁跳变。这个办法不算理论严谨工程上却很管用。6.3 二维 OTSU 的工程化落地建议实际项目里我不会每次都对整张原图跑二维 OTSU。推荐的做法是先做一次基于 ROI 的裁剪只在目标可能出现的区域计算。目标占比小的时候二维直方图里背景类概率碾压目标类类间方差会被背景主导阈值容易往背景一侧偏。另外一个经验是很多场景下可以把二维 OTSU 当作“粗分割”步骤先拿到目标区域再配合形态学开闭运算清理边缘。我处理工业零件尺寸测量时就是先跑二维 OTSU 定位区域再用 Canny 提边缘整体精度比直接用一维 OTSU 高不少。还有一点要提醒二维 OTSU 的计算结果非常依赖图像质量光照变化剧烈的场景建议先做光照校正或归一化。我试过把二维 OTSU 直接用在强阴影图像上效果反而不如一维 OTSU因为阴影区域和目标的邻域均值分布很接近。最后再分享一个调试技巧。如果你不想每次运行都打印阈值可以在程序里把二维直方图保存成图片或 CSV 文件看一眼直方图热力图确认 C、D 区域是不是真的可以忽略。我之前有一批零件图像分割不准就是靠看热力图发现噪声区域占比过高后来调整了窗口大小和 ROI问题很快解决。二维 OTSU 本身不复杂复杂的是在具体场景里把参数和环境调配合适多画图、多对比比单纯调代码有效得多。