企业官网建设流程全解析

简介：这份压缩包提供一套用 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算对，minlength=256这类参数漏掉，空灰度级会被压缩，后面索引对不上，阈值就会整体偏移。

提示：输入必须是单通道 uint8。彩色图直接喂进按单通道写的 Otsu 函数，轻则维度报错，重则把三个通道的灰度混在一起统计，阈值完全失去意义。

2.2 手写版本：一趟循环把阈值搜出来

先看最朴素的写法，目的是把公式和代码一一对应，不追求性能。下面这段可以直接放进任何 Python 环境，只依赖 numpy：

import numpy as np def otsu_naive(gray: np.ndarray) -> int: """朴素版 Otsu：逐阈值计算类间方差，返回最佳阈值。""" # 统计 0~255 每个灰度级的像素数，minlength 防止空级被压缩 hist = np.bincount(gray.ravel(), minlength=256).astype(np.float64) total = gray.size levels = np.arange(256, dtype=np.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。把前缀和提前算好，整个阈值搜索可以压成几次cumsum：

import numpy as np def otsu_vectorized(gray: np.ndarray) -> int: """向量化 Otsu：用前缀累积一次性算出所有阈值的类间方差。""" hist = np.bincount(gray.ravel(), minlength=256).astype(np.float64) total = hist.sum() levels = np.arange(256, dtype=np.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, out=np.zeros(256), where=cum_cnt > 0) # 前景占比与均值：总数减去背景部分 w1 = 1.0 - w0 mu1 = np.divide(mu_total * total - cum_sum, total - cum_cnt, out=np.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/0；sigma_b里可能出现nan，用~np.isfinite统一置 0，防止argmax被污染。这份向量化写法在 2000 万像素的图上通常比朴素版快一个数量级以上，差距主要来自去掉了 Python 层的 256 次循环。

注意：np.argmax遇到并列最大值时返回第一个下标，也就是较小的阈值。如果业务上希望前景多一点，可以改成np.argmax(sigma_b[::-1])再换算，或者对sigma_b做一次平滑再取峰值。这个细节在双峰接近的图上会直接改变输出，不是玄学。

3. 程序怎么接进工程：输入输出、参数与三条调用路径

3.1 文件结构与依赖确认

一份能落地的大津法二值化程序，常见做法是拆成三层：算法核心负责直方图和阈值搜索，图像 IO 负责读图、灰度化、写结果，批处理入口负责目录遍历和并发。核心层只依赖 numpy，IO 层依赖 OpenCV 或 Pillow，批处理层用标准库的argparse和concurrent.futures就够了。下面是这类程序里各模块的职责和输入输出，拿到源码后可以按这张表对号入座：

依赖确认只需要两条命令，先确认 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转灰度后才进 Otsu；GaussianBlur的(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(), minlength=256).astype(np.float64) total = hist.sum() levels = np.arange(256, dtype=np.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, out=np.zeros(256), where=cum_cnt > 0) w1 = 1.0 - w0 mu1 = np.divide(mu_total * total - cum_sum, total - cum_cnt, out=np.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 介于两者之间。下面这张表是我做选型对比时常用的维度：

对照实验的做法是固定同一批图，分别跑三种方法，把阈值、前景占比、人工抽检的误判数记下来。误判数不需要全量标注，每类抽 20 张、每张看 3 个区域就够，重点是看趋势而不是绝对精度。我的血泪经验是：不要用二值图“看起来干净”作为唯一标准，前景占比和阈值稳定性往往更能暴露问题。一张图阈值 128、前景占比 50%，看起来舒服，但换一张同场景图阈值跳到 60，这条流水线就不可靠。

6. 进阶边界：局部 Otsu 与多阈值 Otsu 什么时候该上

全局 Otsu 的天花板在光照不均。一张图左上角亮、右下角暗，全局直方图会被拉成宽峰，单一阈值切下去，亮区前景被吃掉、暗区背景被当成前景。这时候该上局部 Otsu：把图分成若干块，每块独立算阈值，再对块间阈值做双线性插值，避免块边界出现硬接缝。分块尺寸一般取 32 到 64 像素，太小会让平坦块退化成噪声主导，太大又回到全局 Otsu 的老问题。代价是计算量随块数线性增长，2000 万像素的图分 64×64 块，阈值搜索次数是全局版的几千倍，必须走向量化或降采样。

多阈值 Otsu 是另一条路，适合前景本身有多个灰度层次的场景，比如一张图里同时有深色文字、浅色水印和灰色背景。它把类间方差推广到 K 个类，搜索空间从 256 涨到 256^(K−1)，K=3 时还能暴力搜，K 再大就得用动态规划或智能优化算法。我的建议很明确：K=2 用全局 Otsu，K=3 且图像不大可以试多阈值，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和前景占比，最后才敢往流水线里塞。这套习惯帮我省掉了好几次批量翻车的返工。希望帮到你。

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