大津法二值化程序拆解：Otsu原理、工程实现与效果验证

简介这份压缩包提供一套用 MATLAB 编写的大津法二值化程序面向图像处理初学者、课程实验学生以及需要快速验证阈值分割效果的开发者帮助完成灰度图像中前景与背景的自动分离。包内共 1 个文件主体为 .m 脚本压缩包约 1KB代码集中、依赖较少便于直接打开阅读、运行和二次修改。已有 396 人学习适合作为理解经典阈值算法的轻量示例。程序围绕灰度直方图统计、候选阈值遍历、前景与背景平均灰度计算以及类间方差最大化展开读者可据此掌握最佳全局阈值的求解流程并将二值化结果用于文本识别、医学图像分析或工业检测等场景的预处理也能对照公式核查实现细节观察不同图像在噪声和光照变化下的分割表现为后续尝试自适应阈值等方法打下基础。1. 大津法二值化程序为什么“实验证明很好”这句话要拆开看手里这份大津法二值化程序标题写得很直白——“实验证明很好”。我第一反应不是信而是想知道这个“好”是在什么图上、用什么指标、跟谁比出来的。二值化这件事OpenCV 一行cv2.threshold(..., cv2.THRESH_OTSU)就能出结果真正拉开差距的是阈值在直方图双峰不明显时稳不稳、遇到光照不均会不会整片翻白、批量跑几千张时会不会慢到让人放弃。这份程序值得拆的地方是它把 Otsu 的阈值搜索、灰度预处理和输出校验摊开成可读的代码而不是塞进一个黑匣子函数里。适合谁做文档扫描、OCR 预处理、工业缺陷检测这类需要稳定二值化前端的开发者也适合刚接触图像分割、想把类间方差公式和代码对上号的新手。下面按“公式怎么来 → 程序怎么接 → 坑在哪 → 怎么验证”的顺序拆。2. 大津法的数学内核从直方图到类间方差最大化2.1 直方图与类间方差阈值到底在找什么大津法Otsu的核心假设很朴素一幅图里前景和背景的灰度分布大致分成两堆找一个阈值把两堆切开让“堆内尽量像、堆间尽量不像”。数学上不直接算堆内相似度而是等价地最大化类间方差。设灰度级为 0 到 255阈值 t 把像素分成背景类 C0灰度 ≤ t和前景类 C1灰度 t记 ω0、ω1 为两类像素占比μ0、μ1 为两类灰度均值μT 为全图均值类间方差写成σ_b²(t) ω0·ω1·(μ0 − μ1)²遍历 t 从 0 到 255取使 σ_b² 最大的 t 就是 Otsu 阈值。这个公式里有两个容易忽略的细节一是 ω0·ω1 这个权重当某一类像素极少时权重趋近于 0阈值会被人为压向直方图的一侧二是 (μ0 − μ1)² 对两类的均值差敏感如果直方图是单峰或者拖尾很长μ0 和 μ1 的差不会出现明显极大值此时 Otsu 选出的阈值更像是在“矮子里拔将军”。理解这两点后面遇到阈值跳变、整片翻白时就不会慌。实际程序里通常先做灰度化再统计 256 级直方图然后在这张长度为 256 的数组上做遍历。直方图的好处是把逐像素的求和变成按灰度级累加复杂度从 O(W·H·256) 降到 O(W·H 256)。一份能直接跑的大津法二值化程序第一步一定是把bincount或histogram算对minlength256这类参数漏掉空灰度级会被压缩后面索引对不上阈值就会整体偏移。提示输入必须是单通道 uint8。彩色图直接喂进按单通道写的 Otsu 函数轻则维度报错重则把三个通道的灰度混在一起统计阈值完全失去意义。2.2 手写版本一趟循环把阈值搜出来先看最朴素的写法目的是把公式和代码一一对应不追求性能。下面这段可以直接放进任何 Python 环境只依赖 numpyimport numpy as np def otsu_naive(gray: np.ndarray) - int: 朴素版 Otsu逐阈值计算类间方差返回最佳阈值。 # 统计 0~255 每个灰度级的像素数minlength 防止空级被压缩 hist np.bincount(gray.ravel(), minlength256).astype(np.float64) total gray.size levels np.arange(256, dtypenp.float64) best_t, best_var 0, -1.0 for t in range(256): w0 hist[:t 1].sum() # 背景像素数 w1 hist[t 1:].sum() # 前景像素数 if w0 0 or w1 0: continue # 某一类为空方差无意义 mu0 (hist[:t 1] * levels[:t 1]).sum() / w0 mu1 (hist[t 1:] * levels[t 1:]).sum() / w1 # 类间方差w0/total 与 w1/total 是两类占比 var (w0 / total) * (w1 / total) * (mu0 - mu1) ** 2 if var best_var: best_var, best_t var, t return best_t逻辑说明hist是长度为 256 的计数数组下标即灰度级循环里w0、w1是两类像素数除以total得到占比mu0、mu1用加权和除以计数得到均值最后按公式算var。参数说明gray必须是二维np.uint8如果是cv2.imread(..., 0)读进来的图正好满足如果传入 float 图bincount会直接报错需要先astype(np.uint8)并确认值域在 0 到 255。这段代码跑一张 4000×3000 的图循环 256 次、每次对长度 256 的数组切片求和整体在毫秒级但每个阈值都重复求和属于典型的“能跑但不优雅”。朴素版最大的价值是便于调试。当你怀疑阈值不对时在这个循环里打印每个 t 对应的w0、w1、mu0、mu1一眼就能看出是直方图统计错了还是某一类被压得太小。我一般会留一份朴素版作为“参考实现”向量化版本跑出来的阈值和它对不上就说明优化写错了。2.3 向量化与查表把重复计算变成累积和朴素版对每个 t 都重新求前缀和重复计算量是 256×256。把前缀和提前算好整个阈值搜索可以压成几次cumsumimport numpy as np def otsu_vectorized(gray: np.ndarray) - int: 向量化 Otsu用前缀累积一次性算出所有阈值的类间方差。 hist np.bincount(gray.ravel(), minlength256).astype(np.float64) total hist.sum() levels np.arange(256, dtypenp.float64) # 前缀累积cum_cnt[t] 为灰度 t 的像素数 cum_cnt np.cumsum(hist) cum_sum np.cumsum(hist * levels) # 灰度 t 的灰度总和 mu_total cum_sum[-1] / total # 全图均值 # 背景占比与均值除零位置用 where 保护 w0 cum_cnt / total mu0 np.divide(cum_sum, cum_cnt, outnp.zeros(256), wherecum_cnt 0) # 前景占比与均值总数减去背景部分 w1 1.0 - w0 mu1 np.divide(mu_total * total - cum_sum, total - cum_cnt, outnp.zeros(256), where(total - cum_cnt) 0) sigma_b w0 * w1 * (mu0 - mu1) ** 2 sigma_b[~np.isfinite(sigma_b)] 0 # 退化阈值置零 return int(np.argmax(sigma_b))逻辑说明cum_cnt和cum_sum分别对应公式里的背景像素数和背景灰度总和mu0是两者相除前景部分用总量减去背景得到避免再算一次反向累积。sigma_b是一个长度 256 的数组argmax直接返回最佳阈值。参数说明where条件保证某一类为空时不会出现 0/0sigma_b里可能出现nan用~np.isfinite统一置 0防止argmax被污染。这份向量化写法在 2000 万像素的图上通常比朴素版快一个数量级以上差距主要来自去掉了 Python 层的 256 次循环。注意np.argmax遇到并列最大值时返回第一个下标也就是较小的阈值。如果业务上希望前景多一点可以改成np.argmax(sigma_b[::-1])再换算或者对sigma_b做一次平滑再取峰值。这个细节在双峰接近的图上会直接改变输出不是玄学。3. 程序怎么接进工程输入输出、参数与三条调用路径3.1 文件结构与依赖确认一份能落地的大津法二值化程序常见做法是拆成三层算法核心负责直方图和阈值搜索图像 IO 负责读图、灰度化、写结果批处理入口负责目录遍历和并发。核心层只依赖 numpyIO 层依赖 OpenCV 或 Pillow批处理层用标准库的argparse和concurrent.futures就够了。下面是这类程序里各模块的职责和输入输出拿到源码后可以按这张表对号入座模块职责输入输出阈值核心统计直方图、搜索最佳阈值单通道 uint8 数组整数阈值 0~255灰度化BGR/RGB 转灰度、通道校验三通道或单通道图单通道 uint8二值化按阈值生成 0/255 掩码灰度图 阈值二值图批处理入口遍历目录、并发、日志输入目录、参数输出目录 日志依赖确认只需要两条命令先确认 numpy 和 OpenCV 的版本能对上再确认程序入口能被导入python -c import numpy, cv2; print(numpy.__version__, cv2.__version__) python -c from otsu import otsu_vectorized; print(otsu_vectorized.__doc__)第一行打印版本第二行验证核心函数能否被导入。如果第二行报ModuleNotFoundError先确认当前工作目录是不是包的根目录或者用pip install -e .做可编辑安装别急着改代码。参数上唯一需要提前定的是灰度化公式OpenCV 的COLOR_BGR2GRAY用的是 0.299R 0.587G 0.114B如果你手里的图是 RGB 通道顺序直接套这个系数会偏色阈值跟着偏。3.2 命令行跑单张与批量这类程序一般会提供一个命令行入口单张调试和批量跑目录是两种最常用的形态。下面是最小可用的一组调用参数名以你手里的入口为准逻辑是通用的# 单张输出二值图同时打印阈值和直方图峰值位置 python -m otsu.cli --input ./samples/doc.png --output ./out/doc_bin.png --gray-method luma --debug # 批量目录进、目录出保留同名文件4 个进程并发 python -m otsu.cli --input-dir ./samples --output-dir ./out --ext .png,.jpg --jobs 4逻辑说明单张模式走“读图 → 灰度化 → 算阈值 → 二值化 → 写盘”的完整链路--debug会把阈值、直方图双峰位置、前景占比打到标准输出方便判断结果是否合理。批量模式按扩展名过滤输出文件名与输入保持一致--jobs控制进程数。参数说明--gray-method常见取值有luma加权平均、average三通道均值、max最大值法文档类图像用luma彩色背景干扰大的场景可以试max--jobs不要超过物理核数IO 密集时开太多反而拖慢。批量跑之前我一般会先拿 10 张图做小样本试跑确认输出目录结构、命名规则、日志格式都符合预期再放大到全量。原因是批量模式一旦把输出写错目录回滚很麻烦尤其是输入输出目录有重叠的时候。命令行入口如果支持--dry-run先跑一遍只打印不写盘是最省事的后悔药。3.3 Python API 嵌入处理流水线更多时候二值化只是流水线里的一环需要把核心函数嵌进自己的代码。下面这段是典型的接入方式读图、去噪、算阈值、生成掩码、写盘一气呵成import cv2 import numpy as np from otsu import otsu_vectorized def binarize(path: str, out_path: str) - int: img cv2.imread(path, cv2.IMREAD_COLOR) if img is None: raise FileNotFoundError(path) gray cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # 轻度高斯去噪窗口取奇数噪声重时换中值滤波 gray cv2.GaussianBlur(gray, (3, 3), 0) t otsu_vectorized(gray) # 按业务约定决定前景是白还是黑这里前景取白 _, binary cv2.threshold(gray, t, 255, cv2.THRESH_BINARY) cv2.imwrite(out_path, binary) return t逻辑说明cv2.imread默认读成 BGR 三通道cvtColor转灰度后才进 OtsuGaussianBlur的(3, 3)是核大小必须为奇数0表示让 OpenCV 按核大小自动算标准差cv2.threshold里的t就是 Otsu 返回的阈值THRESH_BINARY表示大于阈值的像素置 255。参数说明核大小从 3 起步噪声越重可以加到 5 或 7但核越大边缘越糊OCR 场景一般不超过 5如果图像本身很干净这一步可以省掉省掉的代价是椒盐噪声会把直方图拖平阈值往中间跑。接入流水线时还要注意返回值的用途。binarize返回阈值调用方可以把它记进日志或数据库后续做阈值漂移监控。工业检测里同一批产品图的 Otsu 阈值通常在一个窄区间内波动某天突然整体偏移 20 个灰度级往往意味着光源衰减或镜头脏了这比看二值图本身更早发现问题。4. 常见坑与排查阈值偏移、全白全黑、彩色通道4.1 现象结果图整片白或整片黑这是最高频的翻车现场。现象是输出二值图几乎全白或全黑肉眼看不出任何结构。原因通常有两个一是输入图本身对比度极低比如近乎全黑、全白或大面积纯色直方图退化成一根尖峰Otsu 仍然会返回一个阈值但这个阈值没有区分意义二是彩色图被当成灰度图直接送进bincount三通道展平后灰度值分布被拉宽阈值落在错误位置。解决方式是在算阈值前加一道退化检查def safe_otsu(gray: np.ndarray) - int: # 动态范围过小直接返回固定阈值避免无意义搜索 if int(gray.max()) - int(gray.min()) 10: return int(gray.mean()) return otsu_vectorized(gray)gray.max() - gray.min()小于 10 说明整幅图几乎没有灰度变化这时候任何阈值都切不出结构返回均值至少不会让下游拿到一个极端值。彩色图的解决更简单进 Otsu 之前统一cvtColor并在函数入口用assert gray.ndim 2把维度错误挡在早期。4.2 现象相似两张图阈值跳变很大两张肉眼看起来差不多的图Otsu 阈值一个 80、一个 140二值结果差异明显。原因是 Otsu 对直方图形态敏感当双峰不明显、中间过渡区很宽时类间方差曲线在峰值附近很平微小噪声就能让argmax落到完全不同的位置。这种图用全局 Otsu 本身就不稳解决思路分两层先做直方图平滑或图像滤波把过渡区收窄如果业务允许改用局部 Otsu分块统计再插值。滤波这一步最便宜# 中值滤波对椒盐噪声更有效核大小 3 起步 gray cv2.medianBlur(gray, 3) t otsu_vectorized(gray)中值滤波会把孤立的极值像素替换成邻域中值直方图两端的毛刺被削掉双峰更容易分开。代价是细小的前景结构可能被抹掉字符识别场景要谨慎最好在滤波前后各存一张图对比。4.3 现象彩色图上效果完全不对现象是彩色图跑出来的二值图前景背景乱成一片或者程序直接抛维度错误。原因很直接Otsu 的公式建立在单通道灰度直方图上把 BGR 三通道ravel()后送进bincount得到的分布是三通道混合的结果阈值对任何一个通道都不成立。解决方式是在入口处强制灰度化并把通道顺序确认清楚img cv2.imread(path, cv2.IMREAD_COLOR) gray cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # OpenCV 读进来是 BGR assert gray.ndim 2 and gray.dtype np.uint8如果图是从别的库传进来的比如 PIL 读的是 RGB转灰度要用COLOR_RGB2GRAY用错系数会让红色和蓝色区域的灰度值整体偏移阈值跟着偏。一个简单的自检方法是分别打印 R、G、B 三个通道的均值再和灰度图均值对比差距过大就说明通道顺序或系数有问题。4.4 现象大图批处理慢到不可接受现象是单张图跑得挺快一上批量就慢得离谱CPU 跑满但吞吐很低。原因通常不是 Otsu 本身而是批处理层把图像读写和阈值计算串行化了或者每张图都重复做了一次全量直方图统计。解决方式是把 IO 和计算分开用进程池并发同时确认核心函数走的是向量化版本而不是朴素循环。并发数不要拍脑袋按os.cpu_count()的 70% 起步IO 密集的存储可以再降。4.5 现象二值化后前景背景反了现象是同一批图里有的前景是白、有的是黑下游掩码运算结果对不上。原因是 Otsu 只负责选阈值不负责定义哪一边是前景THRESH_BINARY和THRESH_BINARY_INV的选择取决于业务约定而不同来源的图极性可能相反。解决方式是在配置里固定一个约定并在输出前用前景占比做一次校验如果前景占比超过 80% 或低于 5%大概率极性反了自动翻转或直接报警。这个校验比事后人眼看图可靠得多。5. 效果验证怎么证明“实验证明很好”不是玄学5.1 双峰性与类间方差曲线自检“实验证明很好”要能复现第一步是把类间方差曲线画出来。曲线在正确阈值附近应该有一个明显尖峰峰越尖、峰两侧下降越快说明直方图双峰越清晰Otsu 的结果越可信曲线平坦或者多峰说明这张图不适合全局 Otsu。下面这段把曲线和双峰性指标一起算出来import numpy as np def otsu_diagnostics(gray: np.ndarray) - dict: hist np.bincount(gray.ravel(), minlength256).astype(np.float64) total hist.sum() levels np.arange(256, dtypenp.float64) cum_cnt np.cumsum(hist) cum_sum np.cumsum(hist * levels) mu_total cum_sum[-1] / total w0 cum_cnt / total mu0 np.divide(cum_sum, cum_cnt, outnp.zeros(256), wherecum_cnt 0) w1 1.0 - w0 mu1 np.divide(mu_total * total - cum_sum, total - cum_cnt, outnp.zeros(256), where(total - cum_cnt) 0) sigma_b w0 * w1 * (mu0 - mu1) ** 2 sigma_b[~np.isfinite(sigma_b)] 0 t int(np.argmax(sigma_b)) # 双峰性最大类间方差占全图方差的比例越接近 1 越可分 total_var ((levels - mu_total) ** 2 * hist).sum() / total return {threshold: t, peak_ratio: float(sigma_b[t] / total_var) if total_var 0 else 0.0, curve: sigma_b}逻辑说明sigma_b就是类间方差曲线peak_ratio是峰值占全图总方差的比例常见经验是大于 0.7 说明双峰明显、Otsu 可信低于 0.4 就要考虑换方法。参数说明total_var用全图灰度方差做归一化避免不同图像之间因为对比度差异没法比较。把peak_ratio记进日志批量跑的时候按它排序最低的那几张就是最可能需要人工复核的图。提示peak_ratio不是绝对标准不同场景的阈值要自己标定。文档扫描类图像通常能到 0.8 以上光照不均的现场图可能只有 0.3后者不是程序错了是全局 Otsu 本身不适合。5.2 与固定阈值、自适应阈值的对照实验单看一张二值图说“效果好”没有说服力至少要和固定阈值、自适应阈值放在一起比。固定阈值适合光照稳定的产线自适应阈值适合光照渐变的场景Otsu 介于两者之间。下面这张表是我做选型对比时常用的维度方法关键参数适用场景主要短板固定阈值阈值常数光照、批次高度稳定换批次就要重标全局 Otsu无自动双峰明显、背景干净光照不均时整片翻自适应阈值块大小、C 偏移光照渐变、背景纹理多块大小难调易出块状伪影局部 Otsu分块尺寸、重叠大图局部对比度差异大计算量翻倍块边界要插值对照实验的做法是固定同一批图分别跑三种方法把阈值、前景占比、人工抽检的误判数记下来。误判数不需要全量标注每类抽 20 张、每张看 3 个区域就够重点是看趋势而不是绝对精度。我的血泪经验是不要用二值图“看起来干净”作为唯一标准前景占比和阈值稳定性往往更能暴露问题。一张图阈值 128、前景占比 50%看起来舒服但换一张同场景图阈值跳到 60这条流水线就不可靠。6. 进阶边界局部 Otsu 与多阈值 Otsu 什么时候该上全局 Otsu 的天花板在光照不均。一张图左上角亮、右下角暗全局直方图会被拉成宽峰单一阈值切下去亮区前景被吃掉、暗区背景被当成前景。这时候该上局部 Otsu把图分成若干块每块独立算阈值再对块间阈值做双线性插值避免块边界出现硬接缝。分块尺寸一般取 32 到 64 像素太小会让平坦块退化成噪声主导太大又回到全局 Otsu 的老问题。代价是计算量随块数线性增长2000 万像素的图分 64×64 块阈值搜索次数是全局版的几千倍必须走向量化或降采样。多阈值 Otsu 是另一条路适合前景本身有多个灰度层次的场景比如一张图里同时有深色文字、浅色水印和灰色背景。它把类间方差推广到 K 个类搜索空间从 256 涨到 256^(K−1)K3 时还能暴力搜K 再大就得用动态规划或智能优化算法。我的建议很明确K2 用全局 OtsuK3 且图像不大可以试多阈值K 再大或者图很大优先考虑局部 Otsu 加形态学后处理别硬上多阈值。# 局部 Otsu 的最小骨架分块、逐块算阈值、按块写回 def local_otsu(gray: np.ndarray, block: int 64) - np.ndarray: out np.zeros_like(gray) h, w gray.shape for y in range(0, h, block): for x in range(0, w, block): tile gray[y:y block, x:x block] if tile.size 0: continue t otsu_vectorized(tile) out[y:y block, x:x block] np.where(tile t, 255, 0) return out这段骨架没有做块间插值块边界会有明显接缝落地时要在阈值图上做一次双线性插值再逐像素比较。block取 64 是保守值纹理细的图可以降到 32但要做好阈值抖动的心理准备。从那以后我每次接一个新的二值化程序都会先拿三张图跑一遍一张标准文档、一张光照不均的现场图、一张近乎全黑的废图。三张都过再去看peak_ratio和前景占比最后才敢往流水线里塞。这套习惯帮我省掉了好几次批量翻车的返工。希望帮到你。本文还有配套的精品资源点击获取