企业官网建设流程全解析

简介：一个基于Matlab的二维Otsu自动阈值分割源码，面向图像处理学习者、研究者和工程师，用于解决灰度图像前景与背景自动分离问题。算法采用最大类间方差准则，构建二维直方图并遍历候选阈值，找到使类间方差最大的分割点；相比一维Otsu，二维形式利用像素灰度与邻域均值的联合分布，可更好地抑制噪声和光照不均带来的误分割。压缩包内仅含1个.m文件，整体大小2KB，代码以Matlab脚本形式呈现，注释清晰，涵盖二维直方图统计、阈值搜索、类间方差计算及二值化输出等核心环节，便于阅读、调试和二次开发。已有315人学习下载，适合用作图像分割课程设计、目标检测预处理、医学影像分析或算法对比实验的参考实现，也能帮助深入理解Otsu从一维到二维的拓展思路。

1. 二维Otsu凭什么在分割任务里比一维抗干扰

我做文档二值化时被一维Otsu坑过：灰度直方图明明是双峰，阈值也确实最佳，但扫描件底部阴影那块永远被切成前景。后来换成二维Otsu才明白，问题不是阈值公式错了，而是没有把邻域信息放进去。twodimenOtsu.m 的核心思路是给每个像素增加第二个特征——邻域均值，然后在二维平面搜索使类间分离度最大的阈值向量。这个向量由两个值组成，灰度阈值和邻域均值阈值，它们共同决定像素归属，所以对光照渐变和边缘噪声的容忍度比一维算法高出一截。适合用Matlab做图像分割、批处理二值化，或是想理解二维直方图为什么比一维直方图更稳健的开发者。下文从数学公式、完整源码到调参边界，尽量把这份源码拆开讲透。

2. 二维直方图怎么建，类间方差公式为什么有两个均值

2.1 一维Otsu的局限

一维Otsu在所有灰度级上找一个最佳分割阈值 t，使两个类的类间方差 w0w1(μ0-μ1)^2 最大。关键缺陷很直接：它假设分割边界可以用一个灰度级表示。当背景光照本身就存在斜坡或反射不均时，同一物体在图像左上角和右下角的灰度可以差上30个灰度级，全局单一阈值无法同时满足两处。唯一解决办法是让阈值函数也随空间变化，二维Otsu就是其中一种实现。

但二维Otsu并不是在空间上估计每个位置的阈值，而是把“像素自身灰度”和“局部邻域平均灰度”组成二维特征。如果背景是平滑斜坡，背景像素的邻域均值虽然变化，但仍和自身灰度保持接近；前景物体边缘处，自身灰度和邻域均值会出现较大偏离。于是，特征空间上的类别边界比一维直方图更清晰。

2.2 二维直方图及四个区域

定义灰度级数 L=256。统计函数 h(i,j) 表示灰度值为 i 且邻域均值为 j 的像素个数，其中 j 来自以该像素为中心的窗口平均。除以像素总数得到联合概率 p(i,j)。二维直方图是一个 L×L 矩阵，行列分别对应灰度值和邻域均值。阈值向量 (s,t) 将直方图切成四个区域。

需要特别说明的是，经典实现只把 A、B 当作真正参与分割的类别，C、D 在计算类内均值时被剔除。这样设计的目的，等于在分割结果里保留了一个“不确定带”。代价是 A、B 的概率之和小于等于1。这也是很多从一维公式直接迁移的代码出错的地方。

2.3 分离度公式的来龙去脉

总均值在两个维度上分别计算：μ_iT = Σ_i Σ_j i·p(i,j)，μ_jT = Σ_i Σ_j j·p(i,j)。对给定 (s,t)，背景概率 ω_A 等于矩形 A 的概率和，前景概率 ω_B 等于矩形 B 的概率和。背景均值 μ_iA = (Σ_A i·p(i,j)) / ω_A，μ_jA 同理。如果 ω_A=0，则这组阈值没有意义，直接跳过。

类间分离度常用两维迹的形式：

trσ_B^2 = ω_A[(μ_iA-μ_iT)^2 + (μ_jA-μ_jT)^2] + ω_B[(μ_iB-μ_iT)^2 + (μ_jB-μ_jT)^2]

注意这里不能化简成 ω_Aω_B[(μ_iA-μ_iB)^2 + ...] 的形式，因为 ω_A+ω_B 不等于1。化简后的形式只有当两个类覆盖全部概率时才成立。判断一份源码是否严谨，可以先看它是否对 ω_A、ω_B 分别做归一化。本文给出的 twodimenOtsu.m 采用前一个公式。

3. 用累积直方图把 twodimenOtsu.m 写成常数级查找

3.1 完整源码

需要抄作业的直接拿走，下面几小节讲清楚每一段为什么这么写。代码基于 Matlab R2023b 及兼容版本，需要图像处理工具箱中的 rgb2gray、imfilter；如果只想用基础函数，把 imfilter 换成 filter2 时要自己处理边缘策略。

function [best_s, best_t, bin] = twodimenOtsu(img, r) % 二维Otsu自动分割 % 输入： % img - 灰度图或RGB图 % r - 邻域半径，默认1 % 输出： % best_s - 灰度阈值(0~1) % best_t - 邻域均值阈值(0~1) % bin - 逻辑二值图 if nargin < 2, r = 1; end if size(img, 3) == 3, img = rgb2gray(img); end img = im2double(img); img(isnan(img)) = 0; kernel = ones(2*r + 1) / (2*r + 1)^2; meanImg = imfilter(img, kernel, 'symmetric', 'same'); I = round(img * 255); J = round(meanImg * 255); idx = sub2ind([256, 256], I(:) + 1, J(:) + 1); H = accumarray(idx, 1, [256, 256]) / numel(I); [Gx, Gy] = ndgrid(0:255, 0:255); S = cumsum(cumsum(H, 1), 2); SI = cumsum(cumsum(Gx .* H, 1), 2); SJ = cumsum(cumsum(Gy .* H, 1), 2); mu_iT = sum(sum(Gx .* H)); mu_jT = sum(sum(Gy .* H)); best_s = 0; best_t = 0; best_var = -inf; for s = 0:255 for t = 0:255 w0 = rectS(S, 1, 1, s+1, t+1); if w0 < 1e-12, continue; end w1 = rectS(S, s+2, t+2, 256, 256); if w1 < 1e-12, continue; end mu_i0 = rectS(SI, 1, 1, s+1, t+1) / w0; mu_j0 = rectS(SJ, 1, 1, s+1, t+1) / w0; mu_i1 = rectS(SI, s+2, t+2, 256, 256) / w1; mu_j1 = rectS(SJ, s+2, t+2, 256, 256) / w1; sep = w0 * ((mu_i0 - mu_iT).^2 + (mu_j0 - mu_jT).^2) + ... w1 * ((mu_i1 - mu_iT).^2 + (mu_j1 - mu_jT).^2); if sep > best_var best_var = sep; best_s = s; best_t = t; end end end best_s = best_s / 255; best_t = best_t / 255; bin = (img >= best_s) & (meanImg >= best_t); end function v = rectS(S, i0, j0, i1, j1) if i0 > i1 || j0 > j1 v = 0; return; end i0 = max(i0, 1); j0 = max(j0, 1); i1 = min(i1, 256); j1 = min(j1, 256); v = S(i1, j1) - S(max(i0-1, 1), j1) - S(i1, max(j0-1, 1)) + S(max(i0-1, 1), max(j0-1, 1)); end

这段代码可以直接存成 twodimenOtsu.m 使用。下面解释几个关键点。

3.2 邻域均值窗口和二维直方图的下标映射

r 默认是1，窗口3×3；r=2 是5×5。kernel 全部元素为 1/(2r+1)^2，随窗口大小自适应归一化。imfilter 的边界模式用symmetric而不是默认的零填充，这样图像最外圈像素的邻域均值不会被人为压暗，避免分割结果产生黑边。

sub2ind 那行，把两个0~255的灰度映射到1~256的线性索引。accumarray 把所有落在同一格子的像素计数一次完成直方图，然后除 numel(I) 得到联合概率。如果图像本身已经是逻辑型或uint16，im2double 会把它归一化到0~1，再量化回0~255，这一步顺序不能反过来，否则直接 round(img) 对 uint8 会丢失小数。

3.3 rectS 函数为什么需要四个边界钳制

S 是二维累积和，S(i,j) 表示左上角到 (i,j) 的矩形概率总和。要取任意子矩形，标准容斥公式需要四个角。真正容易出错的是边界：当 s=255 时，区域B的下界是 s+2=257，已经超过矩阵维度。rectS 在入口判断 i0>i1，直接返回0，这样主循环可以减少一层判断。i0-1 可能变成0，所以用 max(i0-1,1) 钳制到1，等效于累积矩阵的第0行全为0。这个细节保证了 s、t 在0到255的闭区间都能被安全遍历。

3.4 主循环里的类间方差与均值分离度

整体均值 mu_iT、mu_jT 在循环前计算，避免在65536次迭代中反复求和。对每个候选阈值，用 rectS 从 SI、SJ 取出 i·p、j·p 的区域和，再除以概率 w0 或 w1，得到两类的二维均值。sep 是按迹形式计算的类间分离度，当两个类的均值向量离总体均值都远，且类概率加权后最大时，就是最可信的分割点。

从实现时序看，for s 在外、for t 在内，列优先遍历可以配合Matlab的缓存访问模式。如果换成 parfor 并行外层循环，要把 best_var 的更新做成 reduction 变量，否则容易产生随机抖动。这个优化留给需要的读者验证。

3.5 参数速查表

4. 分割实战：灰度阈值和邻域均值阈值怎么共同决定像素归属

4.1 二值化的条件必须两维同时满足

很多人在拿到 best_s 和 best_t 后，只用 img > best_s 做二值化。这等于只取了二维Otsu的第一维，丢失了邻域约束。正确方式是：

[bs, bt, bw] = twodimenOtsu(gray, 2); figure; imshowpair(gray, bw, 'montage');

bw 的每个点同时比较 img >= best_s 与 meanImg >= best_t，两个都满足才算前景。若目标本身是暗底亮字，则把两个比较方向都反转，不能只反转其中一个，否则会选中 C 或 D 待定区域。

4.2 与一维Otsu在渐变背景上的对比

测试图像可以自己合成，用 cameraman.tif 加上一个从左到右的亮度斜坡，模拟光照不均。脚本如下。

cam = imread('cameraman.tif'); bg = linspace(0.4, 1.0, size(cam,2)); tilt = double(cam)/255 .* bg; t1 = graythresh(tilt); bw1 = imbinarize(tilt, t1); [bs, bt, bw2] = twodimenOtsu(tilt, 1); figure; subplot(1,3,1); imshow(tilt); title('光照渐变图'); subplot(1,3,2); imshow(bw1); title('一维Otsu'); subplot(1,3,3); imshow(bw2); title('二维Otsu');

一维Otsu在倾斜光照下会把左上角偏暗背景误分为前景，因为那里的灰度整体偏低，被当作人物区域。二维Otsu因为背景的邻域均值也同步偏低，仍然落在 A 区域，人物边缘落在 B 区域，输出更稳定。当然，如果光照倾斜幅度继续加大，二维Otsu也会失效，此时把 r 加大到3也只是延缓而非根治。

4.3 噪声场景下的参数选择

二维Otsu最怕的不是相机白噪声，而是椒盐噪声。一个亮点如果周围也有几个亮点，它就会被算作前景。常见做法是先做3×3中值滤波，再调用本函数。下面表格给出常用的参数组合。

medfilt2 来自图像处理工具箱，处理后噪声点通常消失，但图像会轻微模糊。更激进的做法是先用高斯滤波，但高斯滤波会模糊边缘，导致二维Otsu得到的前景轮廓向内收缩，目视上比中值滤波更明显。所以强噪声优先选 medfilt2。

4.4 用分离度比值判断可信度

主循环里把每次的 sep 存到 sep_grid(s+1,t+1)，循环结束后：

sep_sorted = sort(sep_grid(:), 'descend'); ratio = sep_sorted(1) / max(sep_sorted(2), 1e-12);

ratio 越大越好，一般大于1.5可以放心二值化；如果接近1，说明存在多个近乎等效的阈值，分割结果由噪声主导。这种情况即使做二维Otsu，也不如改用局部自适应阈值，或者先做滤波再重新计算。这个指标在很多开源源码里都没暴露，读取一份新代码时可以自己加两行日志打印出来。

5. 从Matlab到工程部署：二维Otsu的量化加速与移植排错

5.1 灰度级降到64级，循环次数变成4096次

把直方图从256×256降到64×64，遍历代价降为原来的1/16。做法：

levels = 64; I = round(img * (levels - 1)); J = round(meanImg * (levels - 1)); % 后续累积和维度都按 levels 计算，s 从 0 到 levels-1

得到阈值后乘以 255/(levels-1) 再映射回原始灰度。因为二维直方图本身就是计数近似，64级已经能保留叠加峰的位置。对于1080p图像，双循环4096次几乎瞬间完成，适合放在视频处理流水线里。继续降到32级会开始损失小目标，尤其当目标和背景的灰度差不足10个灰度级时，峰位可能漂移。

5.2 移植到OpenCV时要注意矩阵行/列顺序

这是最常见的一个坑。Matlab 的 rectS(S, i, j) 中 i 表示行（灰度），j 表示列（邻域均值）。到了 OpenCV 的 Mat，表达式是 S.at (i,j)，行在前，列在后，本质上一致；但如果你用 C 数组并按下标 [j][i] 声明，就会造成转置。判断方法：二维Otsu得到的两个阈值应该一个接近高灰度区，一个接近低灰度区，如果发现两个值一个很大另一个很小，而且二值图出现沿45°方向的撕裂，优先检查累加矩阵的维度方向。

把 best_var 写入日志，当质量比小于1.2时先做中值滤波再重新分割，比盲目调 r 更可控。这个习惯保留下来，二维Otsu在批量图上基本不会翻车。

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