二维Otsu阈值分割：Matlab实现与工程调参全指南

简介一个用Matlab编写的二维Otsu最大类间方差自动阈值分割源代码面向图像处理初学者、科研人员和需要做图像二值化/前景提取的项目开发者。压缩包内共1个m文件大小约2KB代码精简便于直接阅读或在Matlab工程中复用。算法从二维直方图构建入手将每个像素的灰度值及其邻域均值共同映射到二维平面通过遍历灰度级组合计算前景与背景的类间方差最终筛选出最优阈值完成二值化实现中还涵盖灰度级范围定义、像素权重与均值统计等关键细节可帮助理解一维Otsu向二维扩展的完整思路。对处理目标与背景分布不均、或受噪声干扰的灰度图像这份源码能提供一个可运行的基础示例便于自行改造和对比实验。现有315人学习下载适合想深入掌握阈值分割机制、并希望直接运行验证的读者。1. 一维Otsu在光照漂移下失效二维Otsu用邻域均值把目标从渐变背景里拎出来零件在自然光下拍出来的图像背景灰度从一侧向另一侧慢慢漂移前景和背景的灰度直方图叠在一起。一维Otsu最大类间方差法只按像素自身灰度找阈值遇到这种重叠分布阈值会偏到亮背景一侧把半边背景都当成目标。二维Otsu把“像素自身灰度”和“邻域均值”两个量组成二维直方图用最大类间方差准则在二维平面上找一组阈值s, t目标与背景分别落在四象限的对应区域里光照慢变造成的灰度重叠被邻域信息摊开分割结果比一维Otsu稳得多。这篇笔记适合用Matlab做图像分割、缺陷检测或医学影像处理的工程师。二维Otsu需要掌握三个点二维直方图和四象限划分的数学意义、用累计和矩阵把穷举阈值对的时间降下来、边界与灰度级参数的取舍技巧。下面按原理、流程、源码、调参、实战五个层面把这套方案讲完整给出的代码可以直接落到自己的工程里。2. 二维Otsu的数学原理从一维类间方差到二维直方图的四象限划分2.1 一维Otsu的类间方差公式阈值把灰度分布分成两类再比差异设图像灰度范围 0255第 i 级灰度的像素占比为 p_i。给定阈值 k背景类 C0 为 0k目标类 C1 为 k1255。两类的累计概率和均值如下ω0(k) Σ_{i0}^{k} p_i μ0(k) (Σ_{i0}^{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(p2d1)) 可以直观看到一条沿 45 度方向延伸的高密度带。带子越细说明像素邻域关系越一致阈值分割越可靠带子变粗说明目标内部纹理明显二维Otsu会把一部分纹理细节留在非对角区域分割边界上会带零散小块。理解这个可视化形态比死记公式更容易判断自己的图像适不适合二维Otsu。2.3 四象限划分背景和目标在对角区边界噪声落在非对角区取一对阈值 (s,t)按 i≤s、j≤t 把二维直方图分成四块区域0i≤s 且 j≤t背景类概率 ω0区域1is 且 jt目标类概率 ω1区域2i≤s 且 jt灰度低但邻域均值高背景中的孤立亮点区域3is 且 j≤t灰度不低但邻域均值低暗噪声或边界像素二维Otsu计算类间方差时只考虑对角区域0和区域1区域2和区域3不参与分类。这样定位的是“有空间一致性的区域”单点噪声即使自身灰度值撞进目标区间也会因为邻域均值不匹配而被排除。σB²(s,t) ω0·[(μ0i−μTi)² (μ0j−μTj)²] ω1·[(μ1i−μTi)² (μ1j−μTj)²]遍历 s0..254、t0..254令 σB² 最大的一组 (s,t) 就是最优分割阈值。μ0i、μ0j 是背景类在灰度维和均值维上的质心μ1i、μ1j 是目标类的质心。由于同时约束灰度与邻域均值二维Otsu对局部灰度漂移和孤立噪声比一维Otsu更稳。3. 二维Otsu算法流程与边界处理伪代码、邻域卷积模式和参数设定3.1 伪代码主流程从二维直方图到最优阈值对只需五步二维Otsu的完整流程可以拆成五步任何语言实现都按这个骨架走输入灰度图 I0~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两者对直方图的影响不同边界选项填充方式对二维直方图的影响zero外侧补 0边缘像素邻域均值系统性偏低边界像素偏向区域3目标轮廓外扩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 映射到 1grayLevels构造直方图时可以直接当整数索引用。imfilter 的卷积核要对窗口内像素取平均所以 kernel ones(3)/9。same 保证输出尺寸与原图一致。I_q 减 1 再卷积是因为内部尺度是 1-based减成 0-based 后卷积结果不会产生 1 的偏移最后加 1 回到有效索引范围。3.3 邻域尺寸与灰度级数的相互制约窗口越大背景边界越宽邻域窗口 W 决定二维直方图把多少空间信息纳入统计。W3 时每个像素只看周围 8 个像素对划痕、点状缺陷依然敏感是最常用配置W5 时噪声抑制更强但目标边缘会向外扩展约一个像素分割结果偏“胖”W9 适合大片均值特征明显的图像用在细小纹理图像上会把邻近目标的灰度差异一并抹平。灰度级数对类间方差的影响体现在离散化粒度上。L256 时 p 矩阵里大量格子为 0稀疏结构明显L64 时每个格子平均容纳像素数是原来的 16 倍累计和矩阵更平滑但相邻两个阈值之间的类间方差差值变小最大值位置更容易被噪声扰动。经验值是 L 取 128 对多数 8-bit 图像是精度与速度的平衡点不需要为提速一路降到 32除非只是做预处理粗分割。4. 二维Otsu的Matlab源码实现核心函数、调用脚本与关键代码拆解4.1 核心函数代码基于二维累计和的 twodimenOtsufunction [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 % 一维 OtsuMatlab 自带 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 这类目标占比较大的图像上一维阈值大约落在 90100二维的 s 会偏移到 100120t 更接近目标内部灰度分割图中背景里的孤立白色噪点明显减少。原因是噪点自身灰度够亮但邻域均值偏暗落在 (i≤s, jt) 的区域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)。这套加减对应二维离散积分的容斥关系L256 时把原本 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 时函数内部所有计算都在 1L 的量化尺度上完成返回的 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 矩阵会填充 NaNmax 函数返回的位置可能指向错误阈值。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;量化映射后原始灰度 0255 被压到 1L相邻两个量化级之间相差约 4 个原始灰度级。类间方差在最优阈值附近通常平缓s、t 移动 12 个量化级时分割结果往往只有边缘几个像素的差异。但如果目标与背景灰度差本来就只有 1015 个灰度级降采样到 64 会让两个峰在直方图里连成一个必须回退到 128 或 256 级。grayLevels搜索空间相对256级耗时适用场景25665536 个阈值对1×最终确认、目标与背景灰度差小12816384 个阈值对约 1/4日常批量处理默认项644096 个阈值对约 1/16大图快速筛选、预处理粗分割我一般先按默认 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; endparfor 不允许不同迭代写同一矩阵的同一位置先把每行算成局部向量、再整体写回 sigma能避开冲突。需要并行池时用 parpool 启动否则 parfor 按单worker运行速度没有提升。6.3 分割结果不对时的排查顺序先画直方图再查索引最后查边界结果异常按三层排查。第一层用 imagesc(log(p1)) 画二维直方图看目标和背景是否集中在两个对角簇如果直方图弥散一片问题在图像预处理或邻域窗口选择与这份源码无关。第二层检查直方图统计用的是 1-based 还是 0-based 索引Matlab 中 0 不是合法下标Iq(:)1 这类偏移一旦漏掉直方图数据整体错位一行或一列。第三层核对阈值输出单位s_opt、t_opt 从内部 [1,L] 尺度映射回 [0,255] 时是否乘了 (255/(L−1))映射公式漏掉会导致阈值系统性偏小 12 个灰度级与 graythresh 对比时会出现肉眼不易察觉但实际存在的偏差。本文还有配套的精品资源点击获取