企业官网建设流程全解析

简介：一个用Matlab编写的二维Otsu（最大类间方差）自动阈值分割源代码，面向图像处理初学者、科研人员和需要做图像二值化/前景提取的项目开发者。压缩包内共1个m文件，大小约2KB，代码精简，便于直接阅读或在Matlab工程中复用。算法从二维直方图构建入手，将每个像素的灰度值及其邻域均值共同映射到二维平面，通过遍历灰度级组合计算前景与背景的类间方差，最终筛选出最优阈值完成二值化；实现中还涵盖灰度级范围定义、像素权重与均值统计等关键细节，可帮助理解一维Otsu向二维扩展的完整思路。对处理目标与背景分布不均、或受噪声干扰的灰度图像，这份源码能提供一个可运行的基础示例，便于自行改造和对比实验。现有315人学习下载，适合想深入掌握阈值分割机制、并希望直接运行验证的读者。

1. 一维Otsu在光照漂移下失效，二维Otsu用邻域均值把目标从渐变背景里拎出来

零件在自然光下拍出来的图像，背景灰度从一侧向另一侧慢慢漂移，前景和背景的灰度直方图叠在一起。一维Otsu（最大类间方差法）只按像素自身灰度找阈值，遇到这种重叠分布，阈值会偏到亮背景一侧，把半边背景都当成目标。二维Otsu把“像素自身灰度”和“邻域均值”两个量组成二维直方图，用最大类间方差准则在二维平面上找一组阈值（s, t），目标与背景分别落在四象限的对应区域里，光照慢变造成的灰度重叠被邻域信息摊开，分割结果比一维Otsu稳得多。

这篇笔记适合用Matlab做图像分割、缺陷检测或医学影像处理的工程师。二维Otsu需要掌握三个点：二维直方图和四象限划分的数学意义、用累计和矩阵把穷举阈值对的时间降下来、边界与灰度级参数的取舍技巧。下面按原理、流程、源码、调参、实战五个层面把这套方案讲完整，给出的代码可以直接落到自己的工程里。

2. 二维Otsu的数学原理：从一维类间方差到二维直方图的四象限划分

2.1 一维Otsu的类间方差公式：阈值把灰度分布分成两类再比差异

设图像灰度范围 0～255，第 i 级灰度的像素占比为 p_i。给定阈值 k，背景类 C0 为 0～k，目标类 C1 为 k+1～255。两类的累计概率和均值如下：

ω0(k) = Σ_{i=0}^{k} p_i μ0(k) = (Σ_{i=0}^{k} i·p_i) / ω0(k) ω1(k) = 1 − ω0(k) μ1(k) = (μT − ω0(k)·μ0(k)) / ω1(k)

其中 μT 是整幅图像的灰度均值。类间方差定义为 σB²(k) = ω0·(μ0−μT)² + ω1·(μ1−μT)²，这个值的含义是两类的中心离全局中心越远，说明两类之间拉得越开。遍历 k = 0..254，取 σB² 最大对应的 k 即为最优阈值。

一维Otsu在灰度直方图呈明显双峰时又快又稳，但它只观察像素自身的灰度值，对纹理、孤立噪声和渐变背景没有抵抗能力。灰度均值漂移时，直方图的两个峰宽度加大，波谷位置整体移动，全局阈值就会随之偏向亮背景一侧。这也是二维Otsu被提出的直接动机：在灰度之外引入空间邻域信息，让阈值判据有上下文。

2.2 二维直方图：每个像素用“自身灰度-邻域均值”坐标对建联合分布

二维Otsu对每个像素取两个特征：自身的灰度值 f(x,y)，以及周围一个邻域窗口（通常 3×3）内像素的灰度均值 g(x,y)。把 (f,g) 当作一个二维坐标点，统计整幅图像落在每个格子 (i,j) 的像素个数，再除以总像素数，得到 L×L 的联合概率矩阵 p(i,j)，i 是灰度维度，j 是邻域均值维度。

I = imread('cameraman.tif'); Iq = double(I); Imean = imfilter(Iq, ones(3)/9, 'symmetric'); p2d = accumarray([round(Iq(:))+1, round(Imean(:))+1], 1, [256 256]); p2d = p2d / numel(Iq);

用 imagesc(log(p2d+1)) 可以直观看到一条沿 45 度方向延伸的高密度带。带子越细，说明像素邻域关系越一致，阈值分割越可靠；带子变粗，说明目标内部纹理明显，二维Otsu会把一部分纹理细节留在非对角区域，分割边界上会带零散小块。理解这个可视化形态，比死记公式更容易判断自己的图像适不适合二维Otsu。

2.3 四象限划分：背景和目标在对角区，边界噪声落在非对角区

取一对阈值 (s,t)，按 i≤s、j≤t 把二维直方图分成四块：

• 区域0：i≤s 且 j≤t，背景类，概率 ω0
• 区域1：i>s 且 j>t，目标类，概率 ω1
• 区域2：i≤s 且 j>t，灰度低但邻域均值高，背景中的孤立亮点
• 区域3：i>s 且 j≤t，灰度不低但邻域均值低，暗噪声或边界像素

二维Otsu计算类间方差时只考虑对角区域0和区域1，区域2和区域3不参与分类。这样定位的是“有空间一致性的区域”，单点噪声即使自身灰度值撞进目标区间，也会因为邻域均值不匹配而被排除。

σB²(s,t) = ω0·[(μ0i−μTi)² + (μ0j−μTj)²] + ω1·[(μ1i−μTi)² + (μ1j−μTj)²]

遍历 s=0..254、t=0..254，令 σB² 最大的一组 (s,t) 就是最优分割阈值。μ0i、μ0j 是背景类在灰度维和均值维上的质心，μ1i、μ1j 是目标类的质心。由于同时约束灰度与邻域均值，二维Otsu对局部灰度漂移和孤立噪声比一维Otsu更稳。

3. 二维Otsu算法流程与边界处理：伪代码、邻域卷积模式和参数设定

3.1 伪代码主流程：从二维直方图到最优阈值对只需五步

二维Otsu的完整流程可以拆成五步，任何语言实现都按这个骨架走：

输入：灰度图 I（0~255），邻域尺寸 W，灰度级数 L 1. 对灰度图按 L 级重新量化，记为 I_q 2. 用 W×W 均值卷积核对 I_q 做卷积，得到邻域均值图 I_mean 3. 遍历每个像素，统计 (I_q(x,y), I_mean(x,y)) 落在各格子的次数， 得到 L×L 联合概率矩阵 p 4. 对 p 做二维累计和，同时累计 i·p(i,j)、j·p(i,j) 两个加权累计矩阵 5. 遍历所有 (s,t)，用累计矩阵 O(1) 计算区域0的概率和均值， 代入类间方差公式，记录最大值对应的 (s,t) 输出：最优阈值对 (s,t)、分割结果

第4步是整个算法性能的关键。如果每次更换 (s,t) 都对区域0重新求和，复杂度是 O(L⁴)，256 级灰度下是约 40 亿次浮点加法，Matlab 里直接跑不动。用二维累计和矩阵把区域求和变成四次查表和三次加减，复杂度降为 O(L²)，256×256 的阈值遍历在普通台式机上能在秒级内完成。这也是区分“能用的二维Otsu源码”和“教学演示版源码”的核心分界线。

3.2 边界填充与卷积模式：replicate、symmetric、zero 三选一

邻域均值在图像边界处缺少像素，必须先做填充再卷积。Matlab 的 imfilter 提供多种边界选项，二维Otsu最常用的是 'replicate' 和 'symmetric'，两者对直方图的影响不同：

grayLevels = 256; I_q = round(double(I) / 255 * (grayLevels - 1)) + 1; kernel = ones(3) / 9; I_mean = imfilter(I_q - 1, kernel, 'symmetric', 'same') + 1;

这里 I_q 映射到 1～grayLevels，构造直方图时可以直接当整数索引用。imfilter 的卷积核要对窗口内像素取平均，所以 kernel = ones(3)/9。'same' 保证输出尺寸与原图一致。I_q 减 1 再卷积是因为内部尺度是 1-based，减成 0-based 后卷积结果不会产生 1 的偏移，最后加 1 回到有效索引范围。

3.3 邻域尺寸与灰度级数的相互制约：窗口越大，背景边界越宽

邻域窗口 W 决定二维直方图把多少空间信息纳入统计。W=3 时每个像素只看周围 8 个像素，对划痕、点状缺陷依然敏感，是最常用配置；W=5 时噪声抑制更强，但目标边缘会向外扩展约一个像素，分割结果偏“胖”；W=9 适合大片均值特征明显的图像，用在细小纹理图像上会把邻近目标的灰度差异一并抹平。

灰度级数对类间方差的影响体现在离散化粒度上。L=256 时 p 矩阵里大量格子为 0，稀疏结构明显；L=64 时每个格子平均容纳像素数是原来的 16 倍，累计和矩阵更平滑，但相邻两个阈值之间的类间方差差值变小，最大值位置更容易被噪声扰动。经验值是 L 取 128 对多数 8-bit 图像是精度与速度的平衡点，不需要为提速一路降到 32，除非只是做预处理粗分割。

4. 二维Otsu的Matlab源码实现：核心函数、调用脚本与关键代码拆解

4.1 核心函数代码：基于二维累计和的 twodimenOtsu

function [s_opt, t_opt, segmented] = twodimenOtsu(I, neighborSize, grayLevels) % TWODIMENOTSU 二维最大类间方差法（二维Otsu）阈值分割 % 输入： % I - uint8 灰度图像，范围 0~255 % neighborSize - 邻域窗口尺寸（奇数），默认 3 % grayLevels - 量化灰度级数，默认 256 % 输出： % s_opt, t_opt - 最优阈值对，单位与原始灰度一致 % segmented - 二值分割图，1 为目标，0 为背景 if nargin < 2, neighborSize = 3; end if nargin < 3, grayLevels = 256; end % 1. 灰度量化到 1~grayLevels Iq = round(double(I) / 255 * (grayLevels - 1)) + 1; kernel = ones(neighborSize) / neighborSize^2; Imean = imfilter(Iq - 1, kernel, 'symmetric', 'same') + 1; % 2. 统计联合概率矩阵 p = zeros(grayLevels, grayLevels); idx = [Iq(:), round(Imean(:))]; for k = 1:size(idx, 1) p(idx(k,1), idx(k,2)) = p(idx(k,1), idx(k,2)) + 1; end p = p / numel(Iq); % 3. 累计和矩阵：概率、灰度加权、邻域均值加权 cp = cumsum(cumsum(p, 1), 2); ci = cumsum(cumsum((repmat((1:grayLevels)', 1, grayLevels)) .* p, 1), 2); cj = cumsum(cumsum((repmat(1:grayLevels, grayLevels, 1)) .* p, 1), 2); muT_i = sum(sum(ci)); muT_j = sum(sum(cj)); % 4. 遍历阈值对，查表计算类间方差 sigma = zeros(grayLevels, grayLevels); for s = 1:grayLevels for t = 1:grayLevels w0 = rectSum(cp, 1, 1, s, t); if w0 < 1e-6 continue; end m0i = rectSum(ci, 1, 1, s, t) / w0; m0j = rectSum(cj, 1, 1, s, t) / w0; w1 = 1 - w0; m1i = (muT_i - w0 * m0i) / w1; m1j = (muT_j - w0 * m0j) / w1; sigma(s, t) = w0 * ((m0i - muT_i)^2 + (m0j - muT_j)^2) + ... w1 * ((m1i - muT_i)^2 + (m1j - muT_j)^2); end end [~, idxMax] = max(sigma(:)); [sIdx, tIdx] = ind2sub([grayLevels, grayLevels], idxMax); s_opt = round((sIdx - 1) / (grayLevels - 1) * 255); t_opt = round((tIdx - 1) / (grayLevels - 1) * 255); % 5. 按阈值对原图分割 segmented = double((Iq >= sIdx) & (round(Imean) >= tIdx)); end function val = rectSum(M, r1, c1, r2, c2) % RECTSUM 用累计和矩阵取矩形区域 [r1..r2, c1..c2] 内元素之和 val = M(r2, c2); if r1 > 1, val = val - M(r1 - 1, c2); end if c1 > 1, val = val - M(r2, c1 - 1); end if r1 > 1 && c1 > 1, val = val + M(r1 - 1, c1 - 1); end end

主函数接口把 neighborSize 和 grayLevels 都暴露成可选参数。对 8-bit 灰度图直接调用[s,t,bw] = twodimenOtsu(I)，内部默认 3×3 邻域和 256 级灰度。返回的 s_opt、t_opt 与输入图像灰度同尺度，方便与 graythresh 的结果直接对比。

4.2 调用脚本：与一维Otsu对比，输出分割图和二维直方图

% demo_otsu2d.m I = imread('cameraman.tif'); if size(I, 3) == 3, I = rgb2gray(I); end % 一维 Otsu：Matlab 自带 graythresh level = graythresh(I); bw1 = imbinarize(I, level); fprintf('一维Otsu阈值: %.2f\n', level * 255); % 二维 Otsu [s, t, bw2] = twodimenOtsu(I, 3, 256); fprintf('二维Otsu阈值: s=%d, t=%d\n', s, t); figure; subplot(1,3,1); imshow(I); title('原图'); subplot(1,3,2); imshow(bw1); title({'一维Otsu', ['阈值=', num2str(round(level*255))]}); subplot(1,3,3); imshow(bw2); title({'二维Otsu', ['s=', num2str(s), ', t=', num2str(t)]});

cameraman 这类目标占比较大的图像上，一维阈值大约落在 90～100，二维的 s 会偏移到 100～120，t 更接近目标内部灰度；分割图中背景里的孤立白色噪点明显减少。原因是噪点自身灰度够亮，但邻域均值偏暗，落在 (i≤s, j>t) 的区域2，没有被列入目标类。

4.3 关键代码拆解：累计和矩阵把区域求和降到 O(1)

主循环里的rectSum(cp, 1, 1, s, t)只做四次矩阵索引和三次加减，背后是二维前缀和原理。cp(i,j) 保存矩形 [1..i, 1..j] 内所有概率之和，那么任意矩形 [r1..r2, c1..c2] 的和等于 cp(r2,c2) − cp(r1−1,c2) − cp(r2,c1−1) + cp(r1−1,c1−1)。这套加减对应二维离散积分的容斥关系，L=256 时把原本 6.5 万次 × 6.5 万次的区域累加压缩到常数次查表。

ci 和 cj 的构造同理：repmat((1:grayLevels)', 1, grayLevels) 生成第 i 行全为 i 的矩阵，乘 p 后做累计和，ci(i,j) 就是区域 [1..i,1..j] 内所有“灰度加权概率”之和。区域0的均值向量直接用 ci/cp、cj/cp 计算；区域1的均值则用全局均值减掉区域0的加权贡献再除以 ω1，不需要第二次矩形求和。

grayLevels 不等于 256 时，函数内部所有计算都在 1～L 的量化尺度上完成，返回的 s_opt、t_opt 通过 (idx−1)/(L−1)×255 映射回原始尺度。这样调用方始终面对同样的数值语义：不管内部用什么灰度级穷举，输出阈值都可以直接与 graythresh 的结果做差。分割图 segmented 直接用量化后的 Iq 与 round(Imean) 比较，避免把浮点均值图插值回原尺度。

注意：主循环里if w0 < 1e-6, continue; end不是性能优化，是防止 w1 = 1 − w0 在 w0 接近 1 时出现除零。去掉这个判断，目标占比极小的图像里 sigma 矩阵会填充 NaN，max 函数返回的位置可能指向错误阈值。

5. 运行速度与调参：灰度级降采样、邻域尺寸与二维Otsu的耗时规律

5.1 256级灰度下的时间开销：遍历 65536 个阈值对，查表是收益来源

二维Otsu默认在 256×256 个 (s,t) 上穷举。每个 (s,t) 需要 6 次 rectSum 查表（w0、m0i、m0j 各两次）和约 20 次浮点运算，总计算量在百万次量级，Matlab 里用双层 for 循环跑完整轮通常在 1 秒以内。真正占时间的是卷积计算和直方图统计的像素循环，这部分无法用累计和加速。1080p 分辨率约有 200 万像素，直方图统计循环在纯 M 代码里是大头。

不要用“先遍历灰度再遍历邻域均值”的朴素双循环版本。那种写法在每次更新 (s,t) 时重新累加矩形区域，复杂度 O(L⁴)，256 级时运行时间会到分钟级甚至更久。拿到任何一份二维Otsu源码，先看主循环里有没有调用累计和矩阵，这一条决定代码能不能用于实际图像。

5.2 灰度级降采样：降到64级还是128级，阈值偏移可以接受

把 grayLevels 从 256 降到 128 或 64 是最直观的提速手段：

L = 64; Iq = round(double(I) / 255 * (L - 1)) + 1;

量化映射后原始灰度 0～255 被压到 1～L，相邻两个量化级之间相差约 4 个原始灰度级。类间方差在最优阈值附近通常平缓，s、t 移动 1～2 个量化级时，分割结果往往只有边缘几个像素的差异。但如果目标与背景灰度差本来就只有 10～15 个灰度级，降采样到 64 会让两个峰在直方图里连成一个，必须回退到 128 或 256 级。

我一般先按默认 256 级跑一次拿到参考阈值，再降到 128 级验证阈值变化，若偏移不超过 2 个原始灰度级，就长期使用 128 级。这个验证步骤比盲目选 L 可靠得多。

5.3 类间方差矩阵的形态：最大值附近像一个平台而不是尖峰

把 sigma 矩阵用 surf 画出来，看到的不是尖锐的全局最高点，而是一块弧面平台。灰度级离散化后，目标与背景的边界像素散落在相邻几个格子中，s 和 t 同时增减一两个量化级，区域0和区域1的组成变化很小，类间方差只在小数点后第四位附近浮动。

surf(sigma); xlabel('邻域均值阈值 t'); ylabel('灰度阈值 s'); zlabel('类间方差');

处理这种平台区时，可以适度倾向“低阈值”：在 sigma 值下降不超过 0.5% 的候选点里，取 s、t 更小的一组。低阈值保留更多目标细节，适合缺陷检测这类宁可多检、不可漏检的任务；高阈值则适合希望直接分割出连通区域的场景。两种选择没有绝对优劣，关键是知道 sigma 曲面有这个特征，不要把几个像素级的阈值变化当成算法不稳定。

6. 二维Otsu的三个落地技巧：低对比度增强、parfor并行与失败排查顺序

6.1 低对比度图像先做直方图均衡化，再看阈值是否向期望方向移动

二维Otsu不关心直方图形状，但对比度过低时，联合概率矩阵里大量像素集中在相邻几个格子，类间方差峰值不突出。对图像先执行Ie = histeq(I);再调用函数，灰度范围展开，直方图稀疏化后阈值对更稳定。代价是增强后的灰度不再是物理亮度，阈值不能直接用于原始图像的定量分析。只在“只要分割结果、不要阈值语义”时做这一步。

6.2 parfor并行遍历阈值对：把外层 s 分给多个worker

主循环两层 for 没有数据依赖，内层只依赖累计和矩阵与全局均值，可以并行。把外层换成 parfor 时，注意 sigma 按行写回：

parfor s = 1:grayLevels sigma_row = zeros(1, grayLevels); for t = 1:grayLevels w0 = rectSum(cp, 1, 1, s, t); if w0 < 1e-6, sigma_row(t) = 0; continue; end m0i = rectSum(ci, 1, 1, s, t) / w0; m0j = rectSum(cj, 1, 1, s, t) / w0; w1 = 1 - w0; m1i = (muT_i - w0 * m0i) / w1; m1j = (muT_j - w0 * m0j) / w1; sigma_row(t) = w0*((m0i-muT_i)^2 + (m0j-muT_j)^2) + ... w1*((m1i-muT_i)^2 + (m1j-muT_j)^2); end sigma(s, :) = sigma_row; end

parfor 不允许不同迭代写同一矩阵的同一位置，先把每行算成局部向量、再整体写回 sigma，能避开冲突。需要并行池时用 parpool 启动，否则 parfor 按单worker运行，速度没有提升。

6.3 分割结果不对时的排查顺序：先画直方图，再查索引，最后查边界

结果异常按三层排查。第一层：用 imagesc(log(p+1)) 画二维直方图，看目标和背景是否集中在两个对角簇；如果直方图弥散一片，问题在图像预处理或邻域窗口选择，与这份源码无关。第二层：检查直方图统计用的是 1-based 还是 0-based 索引，Matlab 中 0 不是合法下标，Iq(:)+1 这类偏移一旦漏掉，直方图数据整体错位一行或一列。第三层：核对阈值输出单位，s_opt、t_opt 从内部 [1,L] 尺度映射回 [0,255] 时是否乘了 (255/(L−1))，映射公式漏掉会导致阈值系统性偏小 1～2 个灰度级，与 graythresh 对比时会出现肉眼不易察觉但实际存在的偏差。

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