Otsu大津法二值化程序拆解与分块自适应改造

简介：这份大津法二值化程序面向图像处理初学者、MATLAB 实验者以及需要快速验证全局阈值分割效果的学习者，围绕如何利用类间方差最大化自动寻找最佳阈值、将灰度图转为黑白二值图展开，可服务于文本识别、医学图像分析和工业检测等基础实验场景。压缩包为 rar 格式，共 1 个文件，主体是 1 个 m 脚本，包体约 1KB，体积轻便，便于直接放入 MATLAB 环境阅读与调试。已有 396 人学习，说明其在入门级图像分割练习中具有一定参考价值。通过阅读代码可梳理灰度直方图统计、阈值遍历、前景与背景均值及类间方差计算等关键环节，并对照不同图像观察二值化效果，理解大津法的适用条件与光照、噪声带来的局限，为后续学习自适应阈值等方法打下基础。

1. 拍歪的文档照片和一行阈值：这份大津法程序解决什么问题

上周帮朋友处理一批手机拍的表格照片，要转成黑白稿去做 OCR。我先用固定阈值试了三版，换一张图就大面积失效：亮的那张整片糊成黑块，暗的那张字全被吃掉。换成这套大津法二值化程序之后，一条命令跑完整个目录，肉眼能用的比例从三成涨到八成多。它的定位很朴素：把灰度图自动切成前景和背景两类，阈值不用你手填。适合文档扫描件、表单票据预处理、工业检测里“先把背景分干净再谈别的”这类场景。反过来，光照极不均匀、或者前景只占几个像素的图，全局大津法一定翻车，得上分块或者局部阈值。下面按我拆这份源码的顺序讲：目录怎么分、核心算法怎么落地、怎么跑、哪里会出事、最后怎么改成自己项目要的形状。

2. 工程骨架：从目录结构看一套大津法二值化程序怎么拆模块

拿到一份算法类源码，我先看的不是算法，是目录。目录分得清不清楚，直接决定这东西是“能拿来改”还是“只能拿来跑”。这份程序的骨架属于前者：算法、IO、命令行三层被拆开了，算法层不碰文件系统，所以单元测试可以直接喂数组进去，不用先造一堆图片。

2.1 目录长什么样：算法、IO、命令行三层

otsu-binarize/ ├── main.py # 入口，只做参数解析和调度 ├── core/ │ ├── __init__.py # 对外暴露 otsu_threshold / binarize │ ├── grayscale.py # 通道处理，统一收敛到 uint8 单通道 │ ├── otsu.py # 大津法核心：直方图 + 类间方差 │ └── postprocess.py # 反色、形态学、边缘清理 ├── utils/ │ ├── io_utils.py # 路径扫描、中文路径、批量落盘 │ └── viz.py # 直方图叠加阈值线，调试用 ├── tests/ │ ├── test_otsu.py # 与 OpenCV 内置 OTSU 对拍 │ └── fixtures/ # 合成双峰图，不依赖真实照片 └── requirements.txt

这套分层最关键的一条是core/里没有任何cv2.imread。otsu.py只接受np.ndarray，返回一个整数阈值；binarize再拿这个阈值去做比较。好处很实在：调算法的时候不用管文件名、编码、后缀，跑测试几百毫秒就出结果。我见过不少同学把imread和阈值计算写在一个函数里，最后想验证公式对不对，得先准备十张图。

postprocess.py单独拎出来也是有意的。二值化之后的开闭运算、去小连通域、反色，本质都是“对二值图做修补”，跟阈值怎么算没关系。放一起会导致你想换个阈值算法时，得把修补逻辑一起搬。

2.2 技术栈取舍：OpenCV 打底 + NumPy 手写核心

常被问的一个问题是：OpenCV 一行cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)就能出结果，为什么还要手写一遍？我的答案有两个。第一，手写版能改。比如限定阈值搜索区间、给直方图加权重、把类间方差换成类内方差，这些 OpenCV 都不给你开口子。第二，手写版是黑匣子的对照物。同一张图两边算出来差 1 个灰阶，你就知道问题出在边界处理上，而不是图像本身有问题。

代价也要说清楚。纯 Python 循环遍历 256 个阈值再对每个阈值扫一遍全图，复杂度是 O(256·N)，一张 4000×3000 的图要好几秒。所以核心必须向量化，下面第 3 章会把这部分逐行拆开。

2.3 输入输出约定：批量目录、命名后缀与位深

批量处理的坑，一半出在命名上。这份程序的约定很克制：输入可以是单文件也可以是目录，目录默认不递归；输出统一落到--output指定的目录里，文件名保持原名，后缀换成_bin.png。为什么强制输出 PNG 而不是跟着原后缀走？因为 JPEG 是有损的，二值图存成 JPEG 会在黑白交界处糊出一圈灰边，看起来像算法没做干净，实际是编码器的锅。

from pathlib import Path IMG_EXTS = {".jpg", ".jpeg", ".png", ".bmp", ".tif", ".tiff"} def scan_images(input_path: Path, recursive: bool = False) -> list[Path]: """扫描待处理图片；默认不递归，避免不小心把整个盘扫一遍""" if input_path.is_file(): if input_path.suffix.lower() not in IMG_EXTS: raise ValueError(f"不支持的扩展名: {input_path.suffix}") return [input_path] pattern = "**/*" if recursive else "*" files = [p for p in input_path.glob(pattern) if p.suffix.lower() in IMG_EXTS] return sorted(files) # 排序保证多次运行结果顺序一致 def output_path_for(src: Path, out_dir: Path, suffix: str = "_bin", ext: str = ".png") -> Path: return out_dir / f"{src.stem}{suffix}{ext}"

scan_images里那个sorted不是多余。没有它，glob返回顺序依赖文件系统，批量跑出来的report.csv每次行序都不一样，做前后对比时会误以为结果变了。

位深这一层要在入口卡死。大津法的公式建立在“灰度级可枚举”之上，8 位是 256 级，直接套公式没问题；16 位图有 65536 级，np.bincount(..., minlength=256)会直接抛越界，或者更糟——某些实现里悄悄截断成低 8 位，结果看起来“能跑”，实际上信息已经被砍掉了。常见做法是：16 位图先线性缩放到 8 位再算阈值，或者把直方图的桶数提到 65536，公式本身不用改。

3. 核心算法拆解：直方图累计到类间方差最大化的完整链路

大津法的推导在很多教材里写得很长，落到代码上其实只有四行：归一化直方图、累计概率、累计均值、最大化类间方差。关键在于理解每一步的物理意义，否则调参的时候你不知道该动哪里。

3.1 灰度化与直方图统计：256 个桶怎么来的

输入是彩色图时，第一步是转灰度。OpenCV 的COLOR_BGR2GRAY用的是加权平均，绿通道权重最高、蓝通道最低，这跟人眼对绿光更敏感是一致的。有的人图省事用img.mean(axis=2)做三通道平均，结果红色文字和蓝色背景的亮度差被压平，直方图两个峰贴在一起，阈值就飘了。

import cv2 import numpy as np def to_gray(img: np.ndarray) -> np.ndarray: """统一收敛到 uint8 单通道；带 alpha 的 PNG 先丢掉透明通道""" if img.ndim == 2: return img if img.shape[2] == 4: img = cv2.cvtColor(img, cv2.COLOR_BGRA2BGR) return cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) def gray_hist(gray: np.ndarray) -> np.ndarray: if gray.dtype != np.uint8: raise TypeError(f"直方图只接受 uint8，当前是 {gray.dtype}") # reshape(-1) 把二维拉平，bincount 统计每个灰度级出现次数 return np.bincount(gray.reshape(-1), minlength=256).astype(np.int64)

bincount的minlength=256必须显式给。纯黑图或者纯白图，实际出现的灰度级可能只有一两个，不给minlength就返回长度 1 或 2 的数组，后面按下标访问prob[200]直接 IndexError。这种错误在测试集里很难复现，因为真实照片几乎不会只有一种灰度——但批量处理时，一张全黑的废片就能把整个任务打断。

3.2 类间方差向量化实现：一次 cumsum 走完全部阈值

设阈值为 t，把 0~t 的像素记为背景、t+1~255 记为前景。记 w0(t) 为背景像素占比、μ0(t) 为背景均值、μ1(t) 为前景均值、μT 为全图均值，类间方差是：

σ²(t) = w0·w1·(μ0 − μ1)²

直接按这个式子算，每个 t 都要重新求两次均值，得套两层循环。把它改写一下，用累计量表示：

σ²(t) = (μT·w0(t) − μ(t))² / (w0(t)·(1 − w0(t)))

其中 μ(t) = Σ i·p(i)，i 从 0 到 t，就是累计一阶矩。这样一次cumsum就能把全部 256 个候选阈值的 σ² 都算出来，取最大值对应的下标即可。

def otsu_threshold(gray: np.ndarray) -> int: """大津法：返回使类间方差最大的阈值 t（0~255）""" if gray.ndim != 2: raise ValueError("需要单通道灰度图") if gray.dtype != np.uint8: raise TypeError(f"需要 uint8，当前是 {gray.dtype}") hist = np.bincount(gray.reshape(-1), minlength=256).astype(np.float64) total = hist.sum() if total == 0: return 0 # 空图兜底 prob = hist / total # p(i)：每个灰度级的概率 levels = np.arange(256, dtype=np.float64) omega = np.cumsum(prob) # w0(t) mu = np.cumsum(prob * levels) # μ(t) mu_t = mu[-1] # 全图均值 μT denom = omega * (1.0 - omega) # w0·w1 # 全黑/全白图会让 denom 出现 0，用极小值兜住，避免除零告警 denom = np.where(denom < 1e-12, 1e-12, denom) sigma_b = (mu_t * omega - mu) ** 2 / denom # t=255 时前景为空，没有物理意义，只在 0~254 里找 return int(np.argmax(sigma_b[:255]))

几个参数和边界值得单独说。

sigma_b[:255]这个切片是血泪经验。不切的话，omega[255]恒等于 1，denom被兜底成 1e-12，分子是(μT·1 − μT)² = 0，理论上是 0，但浮点误差可能算出个 1e-13 之类的小量，再除以 1e-12 就放大成 0.1，某些图上真的会赢过正常阈值。切掉最后一项，从根上避免这个玄学。

denom的兜底用np.where而不是denom[denom == 0] = 1e-12。后者的写法在denom是一维数组时没问题，但如果有人把这段代码改成按行处理（比如分块版本），布尔索引会把形状压平，赋值就对不上了。np.where保持形状不变，复制到分块场景更省心。

返回值是一个整数 t，语义是“灰度 ≤ t 判为背景”。如果你的场景是白底黑字，gray > t得到的是背景为 0、字为 255，正好；黑底白字就反了，需要--invert或者直接取反。这个方向问题不看一眼输出图是发现不了的——OCR 引擎吃反色图直接返回空串。

3.3 与 cv2.THRESH_OTSU 对拍：差 1 个灰阶也要查清楚

手写版最怕的不是慢，是“看着对，其实差一点”。所以测试里必须有一条对拍用例，用合成双峰图跑几百组，和 OpenCV 的结果比对。

import cv2 import numpy as np from core.otsu import otsu_threshold def test_match_opencv(): rng = np.random.default_rng(42) for _ in range(200): # 造双峰分布：前景 + 背景各一个高斯，模拟文档图 fg = rng.normal(185, 16, 4000) bg = rng.normal(55, 16, 6000) gray = np.clip(np.concatenate([fg, bg]), 0, 255).astype(np.uint8) gray = gray.reshape(-1, 1) ours = otsu_threshold(gray) cv_thr, _ = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) assert abs(ours - int(cv_thr)) <= 1, f"手写 {ours} vs OpenCV {cv_thr}"

允许差 1 是有原因的。两峰完全对称时，σ²(t) 在相邻两个 t 上的取值可能相等，argmax取第一个，OpenCV 内部实现取后一个，这不是 bug。但如果差到 2 以上，说明边界处理有分歧，得回去看omega的下标起点——有的实现把 t 定义为“前景起始”，差一位就全盘错。

default_rng(42)固定种子这件事别省。不固定的话，某次随机出的分布两峰贴得太近，测试偶发失败，你重跑一次又过了，然后就再也没人管它。种子固定，失败就是稳定失败，能定位。

4. 跑起来：环境安装、命令行参数与批量处理流程

算法看明白了，接下来要让它在自己机器上出图。这一章按“装依赖 → 传参数 → 看结果”的顺序走，每一步都给出错时该看哪里。

4.1 环境准备与依赖版本

python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install -r requirements.txt

requirements.txt里我一般这样写：

numpy>=1.23 opencv-python>=4.7 tqdm>=4.65 matplotlib>=3.6 # 只在 --debug-viz 时导入 pytest>=7.0 # 开发环境用

版本下限不是随便写的。numpy1.23 之前，np.bincount在某些平台上对空数组的行为不太一致；opencv-python4.7 之后imdecode对异常的返回更规范。真出问题时，第一件事是pip list | grep -i opencv看版本，别急着怀疑算法。

有一点提醒：opencv-python和opencv-contrib-python同时装会互相覆盖，症状是cv2能 import 但某些函数找不到。虚拟环境里只留一个。

4.2 命令行参数逐项拆解

import argparse def build_parser() -> argparse.ArgumentParser: p = argparse.ArgumentParser(description="大津法批量二值化") p.add_argument("--input", required=True, help="输入文件或目录") p.add_argument("--output", default="./out") p.add_argument("--recursive", action="store_true") p.add_argument("--invert", action="store_true") p.add_argument("--pre-blur", type=int, default=0, help="高斯核尺寸，0 表示不做预处理") p.add_argument("--report", default="report.csv") p.add_argument("--dry-run", action="store_true") args = p.parse_args() # 高斯核必须是正奇数，偶数核 OpenCV 会直接抛异常 if args.pre_blur and (args.pre_blur % 2 == 0 or args.pre_blur < 3): p.error("--pre-blur 必须是 >=3 的奇数") return p def binarize(gray, invert: bool = False): t = otsu_threshold(gray) # 注意方向：灰度 ≤ t 判为背景 out = np.where(gray > t, 255, 0).astype(np.uint8) return (255 - out) if invert else out, t

--dry-run这个开关看着不起眼，实际是排查效率的关键。拿到一批新图，先 dry-run 一遍，把所有阈值打进 csv，扫一眼阈值分布：集中在 120~160 说明图还算正常；一堆 0 或者一堆 254，几乎可以断定这批图光照有问题或者前景太少，直接去查原图，不用等几百张图写完盘才发现不对。

--pre-blur的取值我一般从 3 起步。核太大（比如 15）会把细笔画一起抹掉，尤其是小字号表格里的竖线；核太小几乎没效果。这个参数的调试顺序是：先 0，看阈值和结果；有问题再上 3；还不行就该考虑光照归一化，而不是继续加大核。

4.3 结果验证：直方图叠加阈值线

出图之后别只看一张。我习惯把直方图和阈值线叠在一起看，双峰清不清楚、阈值落在谷底还是坡上，一眼就能判断这次的结果可不可信。

import matplotlib matplotlib.use("Agg") # 无显示环境也能出图 import matplotlib.pyplot as plt import numpy as np def plot_hist_with_threshold(gray: np.ndarray, t: int, save_path: str): hist = np.bincount(gray.reshape(-1), minlength=256) fig, ax = plt.subplots(figsize=(8, 3)) ax.bar(np.arange(256), hist, width=1.0, color="#888888") ax.axvline(t, color="#c0392b", linestyle="--", label=f"threshold={t}") ax.set_xlim(0, 255) ax.legend() fig.tight_layout() fig.savefig(save_path, dpi=120) plt.close(fig) # 批量跑时一定要关，不关会内存泄漏

阈值线落在两个峰之间的谷底，说明这张图是标准的大津法场景，结果基本可信。阈值线贴在某个峰的内侧，说明双峰重叠严重，二值化会把一部分前景吃掉——这时候不是调参能解决的，得回去处理光照或者上分块。阈值线贴在最左或最右，几本可以判定直方图是单峰的，大津法的前提已经不成立。

plt.close(fig)是必须的。批量跑两千张图，每张都留一个 figure 不关，跑到几百张内存就上去了，最后进程被系统杀掉，日志里啥也没有。

5. 避坑与排查：大津法在真实图片上翻车的 5 个现场

算法本身没什么可错的，错都错在“假设不成立”。这一章记的都是我实际踩过的，按现象、原因、解决三段写，遇到对得上的直接抄。

5.1 现象：一半图正常一半全黑，阈值算出来 30 多

原因：手机拍摄或者扫描仪边缘偏光，整张图左亮右暗。全局大津法只有一个阈值，亮区的背景灰阶可能比暗区的前景还高，模型被拉向亮区，暗区整体判成背景。

解决：先做背景估计再阈值化。常见做法是用大核中值滤波或者形态学开运算估一张背景图，然后gray / background * 128做除法归一化，把缓慢变化的光照压掉，剩下的高频部分再交给大津法。

def flatten_illumination(gray: np.ndarray, ksize: int = 51) -> np.ndarray: # 大核中值滤波抹掉文字，剩下的就是背景亮度场 bg = cv2.medianBlur(gray, ksize) bg = np.where(bg == 0, 1, bg).astype(np.float32) # 防除零 norm = gray.astype(np.float32) / bg * 128.0 return np.clip(norm, 0, 255).astype(np.uint8)

ksize要大于文字笔画宽度，通常是笔画的 3~5 倍。这个数太小，背景场里会残留文字，除法之后文字边缘出现光晕；太大，光照变化的细节被抹平，等于没做。

5.2 现象：阈值稳定跑到 0 或 254，输出要么全黑要么全白

原因：直方图单峰。整张图前景占比不到 2%，或者图本身就是一片均匀背景（拍糊了、曝光过度）。大津法最大化类间方差，在单峰分布下最优解会滑到分布的边缘，这不是算法错了，是前提不成立。

解决：先做一次占比检查。算出阈值后统计(gray > t).mean()，前景占比低于 0.5% 或者高于 60% 就打个标记写进 report，人工复核。另一个办法是限定搜索区间，比如只在 30~220 之间找最大值，把贴边的解排除掉。

5.3 现象：批量跑到一半全是空文件，前面几十张却正常

原因：路径里有中文或空格，cv2.imread在某些版本和平台上直接返回None，不抛异常。代码里没判空，cvtColor(None)才炸，错误信息指向灰度化那一行，跟真实原因差了十万八里。

解决：读写全走imdecode+fromfile的组合，绕开路径编码；同时在读入后立刻判空。

def imread_unicode(path: str): data = np.fromfile(path, dtype=np.uint8) # 按字节读，不经过 OpenCV 的路径处理 if data.size == 0: return None return cv2.imdecode(data, cv2.IMREAD_COLOR) # 调用侧一定要判空，别让 None 流到下游 img = imread_unicode(p) if img is None: failed.append(str(p)) continue

判空之后别忘了把失败的文件名收集起来，跑完统一打印。悄悄跳过最坑，你以为处理了 2000 张，实际只有 1800 张，对账的时候才发现。

5.4 现象：16 位深度图跑完输出一片灰，或者直接抛越界

原因：np.bincount(..., minlength=256)遇到 16 位图里大于 255 的像素值，数组长度会扩到几千，后面的levels数组还是 256 长度，广播失败；某些写法里被静默截断，结果就是把高位的细节全丢了。

解决：入口处卡 dtype，非 uint8 的图显式转换。16 位图做线性缩放而不是右移 8 位，右移会丢掉低位的动态范围。

def to_uint8(gray: np.ndarray) -> np.ndarray: if gray.dtype == np.uint8: return gray if gray.dtype == np.uint16: # 用实际最大值归一，比固定 /257 更稳 mx = int(gray.max()) return ((gray.astype(np.float32) / max(mx, 1)) * 255).astype(np.uint8) raise TypeError(f"不支持的 dtype: {gray.dtype}")

5.5 现象：同一张图重跑两次，阈值差了三五级，边缘毛刺明显

原因：扫描线噪声、椒盐点、JPEG 压缩块效应。这些高频噪声会往直方图的中间区域填像素，把谷底抬高，最大值位置就在附近来回跳。

解决：阈值前加一道轻量去噪，中值滤波对椒盐点最有效，高斯滤波对压缩噪声更合适。但顺序不能反——先滤波再算阈值，绝对不能先二值化再滤波，那样毛刺只是被换成圆角，数量一个没少。

if pre_blur and pre_blur >= 3: gray = cv2.medianBlur(gray, pre_blur) # 椒盐噪声优先中值 t = otsu_threshold(gray)

如果加了滤波阈值还是抖，那就是图的动态范围太窄，gray.max() - gray.min()不到 60 的图基本没救，得回去查采集环节，别在算法上耗。

6. 二次开发：把全局阈值改造成分块自适应

真到项目里，全局大津法八成不够用。最省事的改造方向是分块：把图切成 N×N 的小块，每块独立算阈值，得到一张“阈值图”，再逐像素比较。这样光照不均的图能被局部消化掉。

def tiled_otsu(gray: np.ndarray, tile: int = 64, min_pixels: int = 800): h, w = gray.shape thr_map = np.full((h, w), 128, dtype=np.uint8) for y in range(0, h, tile): for x in range(0, w, tile): patch = gray[y:y + tile, x:x + tile] if patch.size < min_pixels: # 边角小块样本太少，别硬算 continue thr_map[y:y + tile, x:x + tile] = otsu_threshold(patch) # 块与块之间阈值跳变会造成方格状接缝，平滑一下 thr_map = cv2.GaussianBlur(thr_map, (0, 0), sigmaX=max(tile / 4, 1)) return (gray > thr_map).astype(np.uint8) * 255, thr_map

min_pixels那道判断是经验值。64×64 的块只有 4096 个像素，如果里面 90% 都是纯背景，直方图就是单峰，阈值会飘到边上，导致那一小块整片变黑，图上出现一个个黑方块。跳过太小的块，用邻近块平滑后的值兜底，比硬算稳得多。

GaussianBlur的作用是消除接缝。不做这一步，块边界处一边判前景一边判背景，会连成一条纵贯整图的直线，OCR 直接把它当成表格框线。sigmaX取块尺寸的四分之一左右，太大等于退回全局阈值，太小接缝还在。

验证方式我固定用两条。第一条是合成图回归：造一张带人工光晕的渐变背景，中间贴几个已知位置的黑色矩形，跑完检查矩形区域的前景占比是否在 95% 以上，位置是否一个像素都没偏移。第二条是接缝检查：对输出图做一次水平方向的差分，统计整列全变的列数，超过 3 列就说明平滑不够。这两条比肉眼看图靠谱，能进 CI。

一个容易忽略的细节是性能。分块版本对每块都调一次otsu_threshold，每块内部都要跑一遍bincount和cumsum，图大了之后 Python 层的循环开销很显眼。我的做法是先用 1/4 尺寸的缩略图跑一遍，把阈值图放大回原尺寸，再在原图上做比较。阈值本身对分辨率不敏感，缩略图上算出来的 t 和原图上差不了两级，速度却能快十几倍。

从那以后，我每次给图像算法加“看起来更聪明”的改造，都强制先跑一遍合成图回归——光晕渐变、已知位置、固定种子，一样不少。因为在真实照片上，你永远分不清是算法变好了，还是这一批图刚好听话。希望帮到你。

本文还有配套的精品资源，点击获取