1. 为什么一维OTSU不够用从一次分割翻车说起1.1 一维OTSU其实只做了一件简单的事一维OTSU也叫大津法解决的是“给定灰度图能不能自动找到一个阈值把目标从背景里切出来”的问题。它的原理并不复杂把图像像素按灰度值分成两类遍历所有可能的阈值对每个阈值计算两类的类间方差类间方差最大的那个阈值就是最优点。类间方差的计算公式是sigma_b^2 w0 * (mu0 - muT)^2 w1 * (mu1 - muT)^2其中w0、w1是两个类别占全图像的像素比例mu0、mu1是两类像素的灰度均值muT是整幅图像的全局灰度均值。也可以改写成w0 * w1 * (mu0 - mu1)^2这个形式更直观两类均值相差越远或者两类比例越接近五五开类间方差就越大。OpenCV-Python里实现这个算法只需要一行代码ret, binary cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU)这种“遍历所有灰度值找最优分割点”的思路在大部分灰度直方图呈现双峰形态的图像上表现都不错所以它成了经典中的经典。1.2 “只看灰度不看邻居”带来的问题一维OTSU有个致命前提它假设像素的灰度值本身携带了足够的类别信息。可工程实际里这个假设经常站不住。我举一个具体例子图像背景灰度是50目标灰度是200但背景区域有一个孤立的椒盐噪声点灰度是160。单看这个点的灰度值它更接近目标灰度200一维OTSU会毫不犹豫把它归入前景。但是从空间位置看这个点的八个邻居灰度都接近50只有它自己跳变任何有基本图像常识的人都不会认为它是目标。噪声密度一旦上去一维OTSU的二值化结果就会布满孤立小点甚至把真正的目标边界腐蚀得七零八落。做工业缺陷检测或者OCR文本二值化的时候这些噪点会被后续算法当成有效前景连通域分析出来的目标数量翻几倍后处理代码越写越脏。光照不均也是同一个道理同一块物体表面灰度差很大直方图的峰被拉宽甚至重叠一个全局阈值根本没办法同时照顾亮区和暗区这时候OTSU的类间方差最大值对应的阈值往往飘忽不定。1.3 二维OTSU的解决思路给每个像素附加“邻居投票”问题出在特征信息太单薄那解决方向就很明确把空间信息也作为第二个特征喂进判据。二维OTSU的做法是对每个像素计算它邻域窗口内的灰度均值比如3x3窗口然后用二元组当前灰度, 邻域均值来描述这个像素。这样一来背景像素本身灰度低它周围邻居均值也低这个二元组落在二维空间的左下方目标像素灰度高邻居均值也高落在右上方。也就是说正常的目标和背景像素在生成二维直方图时都会聚集在“主对角线”附近。孤立噪声点就不一样了它的灰度可能是160但它的邻域均值接近50这个二元组会掉在“灰度高但邻域均值低”的右下角区域和正常目标相隔很远。我们在这个二维平面上再找一个阈值向量(s, t)就可以把噪声和边缘点单独隔离出去不让它们干扰目标和背景的分类。这就是二维OTSU的基本逻辑不是用更复杂的分类器而是把判据从“一维的点分布”升级成“二维的点分布”让类别之间的分离度计算更鲁棒。2. 二维OTSU原理拆解二维直方图与类间离散度2.1 二维直方图到底长什么样图像灰度值的范围如果取L256每个像素的邻域均值也取0到255的整数那么二维直方图就是一个256x256的矩阵。矩阵元素f(i,j)的含义是全图中同时满足“当前灰度值等于i”且“邻域均值等于j”的像素个数。把它除以像素总数H*W就得到联合概率p(i,j)。p(i, j) f(i, j) / (H * W)这个矩阵的视觉表现很有特点。拿一张光照均匀、噪声不明显的图来说绝大多数像素会集中在主对角线附近因为物体内部灰度均匀某一点的灰度值基本等于它邻域像素的平均值。边缘像素和噪声像素则分布在离主对角线较远的区域但它们的频数通常很小。所以二维直方图本质上就是把“灰度直方图”这张一维统计表升级成了“灰度-局部均值联合统计表”信息量更高但对内存的占用也只是从256个桶变成65536个桶对现代设备完全可以接受。2.2 (s, t)阈值与四象限划分在二维平面上阈值向量(s, t)把整个矩阵切成四个区域A区i s 且 j t灰度低、邻域均值也低主要对应背景B区i s 且 j t灰度高、邻域均值也高主要对应目标C区i s 且 j t灰度低但邻域均值高对应边缘或孤立暗点D区i s 且 j t灰度高的但邻域均值低对应孤立亮点。一维OTSU只有两个类别而二维OTSU实际上把像素分成了四类但C区和D区在理想情况下概率和非常小计算判据时选择忽略它们只保留A区做背景、B区做目标。这个“忽略”的动作非常关键它是抗噪能力的来源孤立的椒盐噪声点落入D区在计算背景和目标统计量的时候直接不参与自然就不会干扰阈值的选取。各类概率和均值向量定义如下。背景概率w0 sum_{i0}^{s-1} sum_{j0}^{t-1} p(i,j)目标概率w1 sum_{is}^{L-1} sum_{jt}^{L-1} p(i,j)背景均值向量u0 (u0i, u0j)目标均值向量u1 (u1i, u1j)全局均值向量uT (uTi, uTj)。这里每个均值都有两个分量一个是灰度分量一个是邻域均值分量。所以二维OTSU处理的不是标量而是一个二维向量。2.3 用迹作为类间离散度的度量一维OTSU的类间方差是标量到了二维空间严格的做法是构造类间离散度矩阵Sb w0 * (u0 - uT)(u0 - uT)^T w1 * (u1 - uT)(u1 - uT)^T这是一个2x2的对称矩阵不能直接比较大小所以取它的迹也就是矩阵对角线元素之和。展开后得到tr(Sb) w0 * [(u0i - uTi)^2 (u0j - uTj)^2] w1 * [(u1i - uTi)^2 (u1j - uTj)^2]这个式子的含义特别直观它同时对灰度维和邻域均值维做了平方误差惩罚两个维度上的类中心偏离整体中心越多迹就越大代表两类在二维空间里分得越开。我们把所有可能的(s,t)都遍历一遍找使tr(Sb)最大的那一对就是最优阈值。相比一维OTSU只优化一个标量方差二维OTSU优化的是一个融合了两个维度信息的标量迹代价是搜索范围从256个候选变成65536个候选。3. Python实现从暴力循环到积分图加速3.1 环境准备与OpenCV-Python依赖实现需要numpy和opencv-python两个核心包。OpenCV-Python的发行包安装名是opencv-pythonimport语句写的是import cv2这个细节经常有人搞混。如果环境里还没有直接pip install opencv-python numpynumpy版本建议1.20以上因为后面会用到np.add.at做二维直方图累加以及np.cumsum做积分图这两个API在新版本里性能更好也更稳定。如果你用的还是老版本numpynp.add.at的行为也没问题只是慢一些。3.2 第一步计算邻域均值图像计算邻域均值最省事的方式是OpenCV的均值滤波cv2.blur。这里有个坑灰度图是uint8类型如果直接对uint8做滤波OpenCV在内部先转成浮点计算但输出类型取决于输入类型结果会四舍五入成uint8丢失精度。所以我会先把图像转成float32再滤波gray cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # 单通道灰度图 blur cv2.blur(gray.astype(float32), (win_size, win_size)) blur np.round(blur).astype(uint8)边界处理也值得说一句。cv2.blur默认borderType是BORDER_REFLECT_101也就是镜像反射填充这比补零更合理。如果用补零图像边缘一圈像素的邻域均值会偏低导致这些本来正常的目标像素被误判成背景。如果你用的函数是平均池化或自定义卷积记得手动指定边界模式。3.3 朴素实现的性能问题拿到二维直方图之后最直接的做法是遍历所有(s,t)对每个候选阈值分别累加A区和B区的概率和均值。每个区域的概率计算如果再用两层循环做累加整体复杂度就是O(L^4)。L256时大约是43亿次内层操作在Python里跑一趟要等接近一分钟这还是在图像本身很小的情况下。我第一次实现时就在这上面吃亏了。实际上二维直方图的区域求和完全可以用二维前缀和积分图来优化。积分图的思想是提前算好从原点(0,0)到每个格子(i,j)的矩形累加和之后任意矩形区域的和只需要查询四个角点再加减即可。同样地不仅要算p(i,j)的积分图还要算ip(i,j)和jp(i,j)的积分图分别用来算灰度分量的均值和邻域均值分量的均值。这样整体复杂度降为O(L^2)也就是65536次查询Python跑起来毫无压力。3.4 积分图加速的完整代码下面给出我这边调试过的完整实现。代码里我用P保存联合概率的积分图用Pi保存灰度加权概率的积分图用Pj保存邻域均值加权概率的积分图。三个积分图都加了一圈padding这样在查询时不需要单独处理边界条件。import cv2 import numpy as np def two_dimensional_otsu(gray, win_size3, L256): 二维OTSU自动阈值分割。 Args: gray: 单通道灰度图uint8 win_size: 邻域均值窗口大小建议3、5、7 L: 灰度级数默认256 Returns: binary: uint8二值图目标为255 best_s: 最优灰度阈值 best_t: 最优邻域均值阈值 blur: 邻域均值图像 h, w gray.shape # 计算邻域均值图像 blur cv2.blur(gray.astype(float32), (win_size, win_size)) blur np.round(blur).astype(np.uint8) # 统计二维直方图 hist2d np.zeros((L, L), dtypenp.float64) np.add.at(hist2d, (gray.ravel(), blur.ravel()), 1) prob hist2d / (h * w) # 构造三个积分图概率、灰度加权概率、邻域均值加权概率 # 所有积分图维度为 (L1) x (L1)P[s][t] 表示矩形 [0,s) x [0,t) 的累加 i_idx np.arange(L, dtypenp.float64).reshape(-1, 1) j_idx np.arange(L, dtypenp.float64).reshape(1, -1) P np.zeros((L 1, L 1), dtypenp.float64) Pi np.zeros((L 1, L 1), dtypenp.float64) Pj np.zeros((L 1, L 1), dtypenp.float64) P[1:, 1:] np.cumsum(np.cumsum(prob, axis0), axis1) Pi[1:, 1:] np.cumsum(np.cumsum(prob * i_idx, axis0), axis1) Pj[1:, 1:] np.cumsum(np.cumsum(prob * j_idx, axis0), axis1) total P[L, L] total_i Pi[L, L] total_j Pj[L, L] # 如果某个区域没有像素直接跳过 if total 0: return np.zeros_like(gray), 0, 0, blur uTi total_i / total uTj total_j / total best_tr -1.0 best_s 0 best_t 0 # 遍历所有可能的阈值对 (s, t) for s in range(1, L): # A区[0,s) x [0,t)概率和直接用P[s][t] for t in range(1, L): w0 P[s, t] if w0 0: continue u0i Pi[s, t] / w0 u0j Pj[s, t] / w0 # B区[s,L) x [t,L)利用积分图的容斥原理 w1 P[L, L] - P[s, L] - P[L, t] P[s, t] if w1 0: continue u1i (Pi[L, L] - Pi[s, L] - Pi[L, t] Pi[s, t]) / w1 u1j (Pj[L, L] - Pj[s, L] - Pj[L, t] Pj[s, t]) / w1 # 计算类间离散度矩阵的迹 tr w0 * ((u0i - uTi) ** 2 (u0j - uTj) ** 2) \ w1 * ((u1i - uTi) ** 2 (u1j - uTj) ** 2) if tr best_tr: best_tr tr best_s s best_t t # 用灰度阈值s直接做二值化 binary np.where(gray best_s, 0, 255).astype(np.uint8) return binary, best_s, best_t, blur if __name__ __main__: # 示例调用 img cv2.imread(test.png) gray cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) binary, s, t, blur_img two_dimensional_otsu(gray, win_size3) cv2.imshow(binary, binary) cv2.waitKey(0)主代码不长只有80行左右。核心步骤就三步算邻域均值图、统计二维直方图、遍历(s,t)找最大迹。后面所有逻辑都围着这三步转。3.5 代码细节说明为什么遍历从s1而不是s0开始因为s0时A区为空w0必定为0没有计算意义。t同理。遍历范围终点取L-1也就是255因为如果s255B区只剩i255这一条线大概率是噪声或少类样本而且它对应的区域最小值也没有意义。实际使用中最优阈值基本不会落在极端边界上所以不需要为了那0.1%的极端情况增加代码复杂度。为什么要用三个积分图而不是直接在循环里累加如果直接在双循环里再去遍历A区和B区复杂度是四次方量级前面已经算过账了。积分图把区域查询变成O(1)操作这是这个实现能在几百毫秒级完成的关键。三个积分图中P只统计概率Pi统计灰度加权概率Pj统计邻域均值加权概率它们共享同一个积分图结构只是累加值不同。代码里用np.add.at来累加二维直方图而不是用双重for循环是因为np.add.at在底层做了索引扁平化速度要快得多。如果你用hist2d[gray, blur] 1numpy的fancy indexing不会执行累加语义只会把最后一个索引对应的位置置1这是另一个容易踩的坑。np.add.at就是为了解决“指定位置的累加”而存在的。4. 实测效果与参数选择4.1 一维与二维OTSU分割结果对比拿一张带高斯噪声的图像测试时差异非常直观。一维OTSU的分割结果里会出现大量独立噪点尤其是在目标与背景灰度差比较大的情况下噪点会被当作前景保留下来二维OTSU因为把邻域均值这个特征拉进来了那些灰度值异常但周围均值正常的孤立点会被划进D区不参与类别统计所以二值图明显干净很多。网上公开的lena图、细胞显微图那些本身质量很高的图像上两种算法差别不大因为像素灰度分布本身就干净。但如果是工业相机拍的带噪声图、手机扫码场景的暗光图二维OTSU的优势会被放大。我印象比较深的一次测试里一张加了高斯噪声的文档图像一维OTSU分割后文字边缘毛刺严重没连接上的笔画碎片很多二维OTSU处理后笔画连贯通顺噪点被控制在极低水平。这个对比说明二维OTSU并不是在所有图像上都大幅优于一维OTSU它针对的是“噪声干扰明显”的场景这一点选用时要心里有数。运行时间上也该有概念。一维OTSU在OpenCV里是C实现的几乎毫秒级返回。二维OTSU的纯Python双重循环版本在256灰度级下需要遍历65536个候选阈值配合积分图查询实测一张1024x1024图像大约需要400~800毫秒取决于机器。如果你的项目对实时性要求高比如30帧处理这个速度就不够看需要考虑缩小灰度级或者用numba/Cython加速。4.2 邻域窗口大小对结果的影响邻域均值窗口的选择直接决定抗噪和保细节的平衡我用表格列出常见窗口的适用场景。窗口大小抗噪能力细节保留适用场景3x3一般最好噪声轻微、目标细小5x5中等较好通用场景推荐先试7x7较强一般噪声明显、目标较大9x9及以上很强差边界膨胀明显大目标、强噪声需谨慎窗口变大时邻域均值会更平滑噪声点更难影响统计结果但目标的细长结构也会被均值模糊掉边缘向外膨胀几个像素是常态。实际项目中我的习惯是先用3x3试如果分割结果还有明显噪点就上5x5还不够再上7x7每次换窗口都要同时检查目标的边缘是否还完整。4.3 阈值对(s,t)与二值化策略很多文章讲完二维OTSU就默认用s对原图做二值化这没错但联合判定也有它的价值。如果你想把C区和D区的杂散点也一并压掉可以用mask (gray best_s) (blur best_t) binary np.where(mask, 255, 0).astype(np.uint8)联合判定的好处是每个孤立噪点要同时满足“灰度低”和“邻域均值低”才可能是背景那些灰度像背景但周围一片亮堂的边缘像素会被划到C区排除。代价是目标边缘细节也可能被丢失。工业场景如果后处理里已经有形态学开闭操作我一般直接用s如果分割结果直接进OCR或检测模型联合判定更稳妥。5. 实际工程中的坑与优化思路5.1 数据精度与归一化实现过程中最容易出问题的地方是数据类型。二维直方图累计的是像素个数一张4000x3000的大图像有1200万像素对应到uint32没有问题但Pi积分图累加的是i*p(i,j)也就是灰度值乘以像素个数最大值可能达到255 * 1200万约30.6亿已经逼近uint32上限42.9亿如果图像更大就直接溢出了。所以稳妥的做法是一开始就用float64做累加虽然内存多占一点但彻底免除溢出风险。邻域均值图像的四舍五入也有讲究。cv2.blur输出的是float32如果不取整直接建二维直方图邻域均值是浮点数没法直接当矩阵下标使用。取整的方式有两种四舍五入和直接截断。四舍五入更接近真实均值直接截断会让均值偏向0我建议用np.round再转uint8。5.2 什么时候二维OTSU也没用二维OTSU不是万金油。目标与背景在灰度上重叠严重的情况它处理不了比如半透明物体叠在相近色背景上无论一维还是二维直方图上两类模式都混在一起判据算不出有意义的峰值。强纹理背景也是同样的问题如果背景本身有很多边缘C区和D区的概率和会明显上升前述“忽略C/D区”的基本假设被破坏分割结果会不稳定。遇到这种图像更合理的方向是先做光照校正、配准或者纹理抑制把图像质量拉回来之后再用OTSU类算法。还有一种情况是光照梯度非常强比如暗角严重的镜头图像左上角和右下角同一物体的灰度差可能超过100。此时一维OTSU的阈值会偏向像素多的区域二维OTSU虽然因为邻域均值同步偏移效果比一维好但依然不是根治方案。工程上我会先做背景估计和相减再套阈值分割。5.3 性能优化缩小灰度级如果摄像头分辨率高、帧率高又必须在Python里跑二维OTSU可以把灰度级从256压到64甚至32。做法很简单取图像灰度值的高6位或高5位即gray 2邻域均值图像也同样处理。二维直方图随之缩小到64x64或32x32遍历次数从65536次降到4096次甚至1024次耗时可以减少一个数量级。灰度级压缩后阈值精度会有轻微损失但对绝大多数分割任务没有可感知的影响。更进一步如果想要更快的运行速度可以用numba的jit装饰器重写双循环部分或者把双循环里计算迹的逻辑用cupy放到GPU上。不过我的经验是灰度级压到64之后纯Python实现就已经够应付大部分离线处理场景真正需要毫秒级响应的场景更推荐用C来实现这个算法。5.4 扩展思路从二维到三维OTSU二维OTSU把灰度维和邻域均值维放进判据那自然也可以再加入第三个特征比如局部梯度幅值值或局部纹理方差形成三维OTSU。三维OTSU的判据从2x2矩阵的迹变成3x3矩阵的迹信息更丰富但候选阈值的搜索空间从256^2暴增到256^3也就是上千万个候选纯遍历不现实通常要用智能优化算法去搜索。这个方向学术文章不少工程落地则要评估收益是否值得。对抗噪要求极高的场景我更推荐先用二维OTSU跑一版结果观察残留噪声的形态特征再有针对性地加形态学处理或后滤波这样往往比直接上三维更可控。最后再分享一个项目里的小经验。我在做元器件外观缺陷检测时打光环境总是不够理想相机采集的图片经常带随机的反光噪声。一开始一维OTSU做出来的二值化图后续连通域分析经常多出上百个小区域。换成二维OTSU之后噪声引起的孤立区域少了一个量级很多针对噪点的后处理规则都可以直接删掉了。代价是代码从一行变成几十行但换来的是稳定。如果你也被这类自动阈值问题困扰建议先把上面这段代码跑通再根据你自己的图像特性去调邻域窗口和判定策略大概率能省下一大堆后处理时间。