简介二维Otsu最大类间方差算法的Matlab实现源码面向正在学习图像分割或需要自动阈值选取的开发者与研究人员。与常见一维Otsu不同二维版本在像素灰度基础上引入邻域平均灰度作为第二维度能更稳健地处理含噪声、前景背景分布不均的图像。压缩包仅包含1个m文件大小约2KB核心代码集中在twodimenOtsu.m中结构紧凑便于直接阅读、调试与二次开发。实现中覆盖了二维直方图构建、灰度级范围设定与阈值初始化、前景/背景权重及均值统计、类间方差公式求解、最优阈值遍历以及最终二值化输出等完整流程并通过Matlab自带的直方图与图像处理函数串联各步骤。目前已有315人浏览学习适合初学者逐段理解算法原理也可作为实验对比或实际项目的轻量参考。通过阅读源码能清晰看到二维Otsu与一维版本的统计区别也方便直接调整阈值参数并套用到自己的图像数据上是快速上手自动阈值分割的一份实用示例。1. 从一维Otsu到二维Otsu为什么单靠像素值不够做图像分割时Otsu是大部分人想到的第一个自动阈值算法。它不需要人工调参遍历所有灰度级找出类间方差最大的那个值把图像分成前景和背景。这个思路在一类图像上表现不错目标和背景灰度差异明显、分布相对集中的简单场景。可一旦图像里出现光照不均、噪声或边缘过渡带一维Otsu的弱点就暴露了——它只考虑每个像素自身的灰度完全没有利用像素与周围邻居的关系。结果往往是阈值偏移浅色目标被割掉一块或者背景里的噪点被当成前景。二维Otsu把单个像素的灰度值和它的邻域均值组成一个二元组在二维直方图上做类间方差最大化。这样既保留了一维Otsu的自动阈值特性又额外考虑了空间上下文信息对噪声和光照变化的鲁棒性明显提升。代价是计算量从256次累加膨胀到65536个格子上的二维遍历但今天的中端PC上跑一张512x512图像毫秒级能完成这个代价是完全可以接受的。本文用MATLAB直接给你一套可复现的二维Otsu实现包含二维直方图构建、类间方差求解、阈值反查三个核心环节并给出参数调节方向和验证方法。无论你是要做细胞图像分割、工业缺陷检测还是只是想把手里的分割算法基线往上提一点这套代码都能直接改着用。2. 二维Otsu算法模型从一维公式到二元组的类间方差推导2.1 一维Otsu的数学回顾与缺陷定位一维Otsu遍历灰度级t把像素分为C0灰度≤t和C1灰度t两类。设总像素数为N灰度级为0到L-1灰度i的出现概率为pi n_i / NC0类出现的概率ω0和均值μ0、C1类的ω1和μ1则有类间方差σ^2(t) ω0 * ω1 * (μ0 - μ1)^2σ^2最大时对应的t即为最佳阈值。这个公式的前提是目标和背景在灰度轴上存在清晰的分界。但实际图像中光照不均匀会让相同目标在不同区域的灰度出现漂移噪声会让个别像素偏离真实灰度值。一维Otsu只对这一个个孤立的灰度值做统计当灰度分布出现多个波峰重叠的时候分割结果会振荡有时甚至会把阴影区域整体误判为前景。我实际调试分割算法时发现一个规律当一维Otsu的分割结果出现目标边缘参差不齐、内部出现孔洞的时候往往不是阈值t本身的问题而是模型本身有缺陷——它缺少一个维度。2.2 二维直方图的定义与三个分区二维Otsu把每个像素用一个二维向量表示x f(i, j) // 像素灰度 y g(i, j) // 像素邻域均值灰度通常取该像素周围3x3或5x5邻域的平均灰度。这样每个像素就映射到二维平面上一个坐标点(x, y)统计所有像素落在每个(x, y)格子上的频数得到L×L的二维直方图记为p(x, y)。在几何分布上二维直方图的点会沿着对角线聚集。原因是自然图像中像素灰度与邻域均值通常高度相关——平坦区域两者接近相等边缘处像素灰度偏离邻域均值则散落在对角线两侧。因此二维直方图天然分成三个区域分区位置含义对角区x ≈ y目标或背景的内部像素灰度均匀边缘区x 与 y 差异大目标边界、噪声点、细纹理交叉区两簇之间的过渡带背景与目标之间的灰度过渡基于这个分布二维Otsu用阈值向量(s, t)把直方图划分为四个象限对角线上两个象限分别对应背景和目标的主体区域反对角线上的两个象限对应边缘和噪声区域。类间方差的计算不仅要让背景和目标各自内部尽量“聚拢”还要让两个分区离得远同时把边缘点对分区归属的干扰降到最低。2.3 二维类间方差公式与求解目标设阈值向量为(s, t)落在背景区域A的概率为ω0 Σ_{x0}^{s} Σ_{y0}^{t} p(x, y)目标区域B的概率为ω1 1 - ω0。背景区域中像素灰度的均值向量μ0 ( Σ_{x0}^{s} Σ_{y0}^{t} x * p(x, y) / ω0 , Σ_{x0}^{s} Σ_{y0}^{t} y * p(x, y) / ω0 )整体灰度均值向量μT (μx, μy)可以在遍历前一次算好。二维Otsu的类间方差定义为离散度矩阵的迹σ^2(s, t) ω0 * [ (μ0x - μx)^2 (μ0y - μy)^2 ] ω1 * [ (μ1x - μx)^2 (μ1y - μy)^2 ]遍历所有(s, t)组合取σ^2最大时对应的(s, t)即为最佳阈值向量。注意这里的均值向量计算用的是概率加权求和不是直接对图像灰度求和这是和一维Otsu在实现上最容易被忽略的区别。用矩阵迹的好处是避免计算协方差矩阵的完整特征值把二维向量的离散度压缩成一个标量。实际使用中特征值分解的排序结果和迹的排序结果在绝大多数图像上几乎一致但迹的计算简单很多对实时性要求高的场景更友好。3. 二维直方图的MATLAB构建与查表法加速3.1 邻域均值灰度图的生成方法二维Otsu的第一步是生成每个像素的邻域均值灰度图这在MATLAB里可以用卷积一步完成function [f, g] compute_gray_pair(img, winSize) % img 输入灰度图像 double 类型范围 0~255 % winSize 邻域窗口大小一般为 3 或 5 f img / 255 * 255; % 保持数值范围不缩放直接用原灰度 kernel ones(winSize) / (winSize^2); g conv2(img, kernel, same); endconv2的same选项保证输出尺寸和原图一致。边界处卷积会引入误差因为图像外部的像素默认补零参与平均导致边缘的邻域均值被拉低。对普通图片来说边缘误差只影响最外侧一圈像素可以忽略。但如果目标是医学影像这类边缘内容非常关键的图我会改用padarray做replicate填充后再卷积代价是稍微多几次运算。3.2 概率矩阵直方图统计的两种方式有了f和g接下来统计二维直方图。最常见也最不容易出错的方式是用三重循环累加function [hist2d, Ps] build_2d_hist(f, g, L) % L 灰度级数例如256 hist2d zeros(L, L); [M, N] size(f); for i 1:M for j 1:N x f(i, j) 1; % MATLAB索引从1开始 y g(i, j) 1; hist2d(x, y) hist2d(x, y) 1; end end Ps hist2d / (M * N); end三重循环的问题在MATLAB里很突出512x512图像意味着262144次迭代即使内部操作很轻也要几十毫秒。换成accumarray可以快一个数量级function [hist2d, Ps] build_2d_hist_fast(f, g, L) f1 f(:) 1; g1 g(:) 1; idx sub2ind([L, L], f1, g1); hist2d accumarray(idx, 1, [L*L, 1]); hist2d reshape(hist2d, L, L); Ps hist2d / numel(f); endsub2ind把二维坐标映射成一维索引accumarray一次性完成累计最后reshape回二维概率矩阵。这里有一个隐藏坑灰度值可能是double浮点比如0.5这种非整数直接用会生成非法索引。所以我在双线性插值或者几何变换之后做Otsu时会先确认输入灰度已经归一化到0~255且取整或者直接用uint8类型输入1转索引时就安全了。3.3 积分图技术在遍历前的预处理价值在遍历256×256个(s, t)组合找最优阈值时如果每步都从头累加ω0和μ0那整体复杂度是O(L^4)256灰度级下就是43亿次运算。MATLAB里跑一次要几分钟完全不可用。解决方式是用积分图。对概率矩阵Ps和加权概率矩阵x·Ps、y·Ps分别建立二维前缀和function S integral_image(M) S cumsum(cumsum(M, 1), 2); end区域(x1:x2, y1:y2)的和通过region_sum S(x2, y2) - S(x1-1, y2) - S(x2, y1-1) S(x1-1, y1-1)在O(1)时间内得到。对每个(s, t)背景区域的三个量——概率和、x加权和、y加权和——分别用一次积分图查询搞定总复杂度降到O(L^2)。L256时65536次查询MATLAB内循环一次就算完了。4. 可运行的二维Otsu MATLAB函数与参数说明4.1 完整实现代码下面给出一个可直接贴进MATLAB运行的函数输入灰度图uint8或double均可输出阈值向量(s, t)和分割二值图。function [s, t, bw] otsu_2d(img, winSize) % img 输入灰度图像uint8或double范围0~255 % winSize 邻域窗口尺寸推荐3或5 % 输出s, t为二维阈值bw为分割二值图 if isinteger(img) img double(img); end L 256; % 1. 生成像素灰度与邻域均值灰度对 kernel ones(winSize) / (winSize^2); g conv2(img, kernel, same); f img; % 2. 构建二维直方图快速accumarray版本 f1 round(f(:)) 1; g1 round(g(:)) 1; idx sub2ind([L, L], f1, g1); hist2d accumarray(idx, 1, [L*L, 1]); hist2d reshape(hist2d, L, L); Ps hist2d / sum(hist2d(:)); % 3. 计算整体均值向量 [xx, yy] meshgrid(0:L-1, 0:L-1); muT_x sum(sum(xx .* Ps)); muT_y sum(sum(yy .* Ps)); % 4. 建立三个积分图 P_int cumsum(cumsum(Ps, 1), 2); XP_int cumsum(cumsum(xx .* Ps, 1), 2); YP_int cumsum(cumsum(yy .* Ps, 1), 2); % 5. 遍历所有阈值组合找最大类间方差 maxSigma 0; s 0; t 0; for ts 1:L % 遍历前景/背景灰度分界 for tt 1:L % 遍历邻域均值分界 w0 P_int(ts, tt); if w0 0 || w0 1 continue; end mu0_x XP_int(ts, tt) / w0; mu0_y YP_int(ts, tt) / w0; mu1_x (muT_x - XP_int(ts, tt)) / (1 - w0); mu1_y (muT_y - YP_int(ts, tt)) / (1 - w0); sigma w0 * ((mu0_x - muT_x)^2 (mu0_y - muT_y)^2) ... (1 - w0) * ((mu1_x - muT_x)^2 (mu1_y - muT_y)^2); if sigma maxSigma maxSigma sigma; s ts - 1; t tt - 1; end end end % 6. 阈值分割灰度大于s且邻域均值大于t则判为目标 bw (f s) (g t); end4.2 分段逻辑说明第1步的卷积核ones(winSize)/(winSize^2)是均值滤波的标准做法。winSize3时每个像素的邻域均值来自周围8个邻居加自身共9个像素的算术平均窗口越大对噪声抑制越强但边缘会变得更模糊。第2步用round()处理灰度值取整防止sub2ind生成非整数下标。accumarray的第一个参数是下标向量第二个参数1表示每次累加1输出长度L×L的列向量reshape回矩阵就是二维直方图。第3步用meshgrid生成两个坐标矩阵xx .* Ps表示“灰度值 × 概率”。这里乘法运算量是L×L65536次MATLAB一次向量化操作完成没必要再用循环。第4步建立三个积分图每个都是256×256矩阵。cumsum的第二个参数1表示按行累积2表示按列累积嵌套调用就是二维前缀和。第5步的双重循环是核心。注意w0取的是P_int(ts, tt)这正好是从(1,1)到(ts,tt)矩形区域的概率和对应二维直方图的背景区域。mu1_x的计算用的是整体均值减去背景区域的加权和再除以前景概率。类间方差公式里(mu0_x - muT_x)^2衡量背景相对整体的偏离程度muT_x - XP_int(ts,tt)则是通过占比推导前景的累计加权和。第6步分割条件用了逻辑与像素灰度本身够高同时它周围的邻域均值也够高才会被判为目标。这个双重条件就是二维Otsu比一维Otsu对孤立噪声点更鲁棒的关键——单点灰度再亮如果周围区域整体偏暗也不会被误判。4.3 复杂度分析与运行时长观察遍历部分虽然是双重循环但每次迭代只做常数次乘法和加法总共65536轮MATLAB R2023b上跑完大约0.3到0.8秒。耗时的主要部分在accumarray和三个积分图的构建这部分是向量化操作几乎可以忽略。如果图像尺寸很大比如4000×4000accumarray构造的索引向量有1600万元素内存占用约128MB。这种规模建议先对图像做降采样再求阈值然后把阈值映射回原图做分割因为Otsu本质是统计量求阈值缩放到原图1/4大小后阈值偏移通常不超过2个灰度级。5. 二维Otsu的关键参数调优与三类常见坑5.1 邻域窗口winSize的选择依据winSize直接影响g矩阵的平滑程度和运算量。给它不同的值分割结果会有肉眼可见的差别。winSize噪声抑制能力边缘保留能力适用场景3中等较好细胞图像、文字提取、轻微噪声5较强中等工业检测、含椒盐噪声的图像7以上强明显变差纹理背景、粗粒度目标分割我一般先用3观察分割结果如果背景噪声破碎点太多再往上加。winSize增大到5以上时边缘被过度平滑细长目标容易断裂5x5是我在实际项目里用得最多的窗口。5.2 灰度范围压缩的加速技巧如果图像实际灰度范围很窄比如工业相机拍到的8位图只有60~180段有像素可以先把灰度压缩到64级或128级再算Otsu。MATLAB实现scale 2; % 128级 img_q floor(img / scale); [s_q, t_q, ~] otsu_2d(img_q, 3); s s_q * scale; t t_q * scale;灰度级从256降到128遍历次数从65536降到16384提速约4倍阈值误差在压缩粒度范围内此处为2个灰度级。对于在线检测这类对延迟敏感的场景代价可接受。提示压缩前先确认图像没有明显的直方图截断。如果min灰度接近0且max接近255说明动态范围本身已经满格压缩会导致信息损失慎用。5.3 一维Otsu与二维Otsu的结果差异核对方法写完函数后先别急着往项目里集成。用同一张图把一维Otsugraythresh和二维Otsu的结果并排显示核对差异是否符合预期% 一维结果 t1 graythresh(uint8(img)) * 255; % 二维结果 [s2, t2, bw2] otsu_2d(img, 3); bw1 img t1; % 统计差异区域占比 diffRatio sum(sum(xor(bw1, bw2))) / numel(img); fprintf(一维阈值: %.1f, 二维阈值: (%d, %d), 差异像素占比: %.2f%%\n, ... t1, s2, t2, diffRatio * 100);正常预期是差异像素占比在5%到20%之间。如果超过30%优先怀疑两件事一是邻域均值g的范围没有和f对齐比如g的值被conv2缩放了二是循环索引偏移弄错了。这里有个容易出现的问题conv2默认补零边界一圈像素的邻域均值小于内部像素。当阈值靠近图像两侧时边界上的像素可能因为均值偏低被判为背景。解决方式是把conv2换成imfilter(img, kernel, replicate)边界用复制扩展而不是补零。5.4 分割结果出现整片反白的排查思路如果bw全为1说明阈值向量(s, t)落在(0,0)附近。可能的原因有三个按出现频率排序第一输入图像是反相的即背景亮目标暗。二维Otsu公式默认前景是灰度高的区域遇到灰度反转图像时(f s) (g t)判断条件要改成(f s) (g t)。第二直方图中背景区域的w0在遍历中只出现过很窄的范围比如w00.99以上导致方差被极端占比主导。检查maxSigma对应位置是否在二维直方图角落如果是说明输入的图像对比度太低或者根本就是单峰分布。第三winSize设置的均值滤波泄漏了灰度边界导致g值整体往中间靠拢阈值向量跟着偏移。排查手段是在MATLAB里画出二维直方图和最佳阈值点的位置figure; imagesc(log(Ps 1e-6)); colormap jet; colorbar; hold on; plot(t 1, s 1, r*, MarkerSize, 15);阈值点应该落在两簇分布交界的谷底。如果红点完全不在谷底附近检查代码里的索引偏移和加权矩阵乘法这两个步骤是错误率最高的地方。6. 二维Otsu的边界情况与自适应改进技巧6.1 不均匀光照下的分块策略对整幅图做一次二维Otsu遇到明显的光照渐变时效果仍然打折扣。改进方向是把图像切块每一块独立做Otsu再用双线性插值把阈值矩阵映射回全图。以256×256的块为例每块大约包含65536个像素足够保证直方图统计的可靠性。块与块之间可能出现阈值跳变导致分割后的二值图出现块状痕迹。解决办法是给每块扩展8像素的重叠区用重叠区的像素一起统计直方图这样相邻块阈值变化会平滑很多。function [sMap, tMap] otsu_2d_blockwise(img, blockSize, winSize) [M, N] size(img); sMap zeros(ceil(M/blockSize), ceil(N/blockSize)); tMap zeros(ceil(M/blockSize), ceil(N/blockSize)); rowBlocks 1:blockSize:M; colBlocks 1:blockSize:N; for i 1:length(rowBlocks) for j 1:length(colBlocks) r1 max(1, rowBlocks(i)); r2 min(M, r1 blockSize - 1); c1 max(1, colBlocks(j)); c2 min(N, c1 blockSize - 1); [sMap(i,j), tMap(i,j)] otsu_2d(img(r1:r2, c1:c2), winSize); end end % 上采样回原图尺寸 s_full imresize(sMap, [M, N], bilinear); t_full imresize(tMap, [M, N], bilinear); bw (img s_full) (imfilter(img, ones(winSize)/(winSize^2), replicate) t_full); end这种分块加上采样的模式比单独做一次全局Otsu在光照渐变的文档图像上分割准确率通常能提升10个百分点以上。代价是计算时间线性增加块数越多越明显。6.2 定量验证分割质量的一个低成本指标确定阈值后除了肉眼看图还可以算一个分割对比度指标来量化结果好坏bgPixels img(~bw); fgPixels img(bw); contrast abs(mean(fgPixels) - mean(bgPixels)) / sqrt(var(double(fgPixels)) var(double(bgPixels)));对比度值大于2说明目标与背景分离良好小于1意味着分割几乎失效。批量测试不同winSize时用这个指标可以快速筛选出最优参数不必一张一张看图。二维Otsu的另一个常见变体是把类间方差公式里的两个维度做加权比如给邻域均值维度的方差乘0.7增强对灰度维度变化的敏感度。这相当于在分割质量和噪声抑制之间调平衡。手头图像噪声水平不高时加权到0.8~0.9能让边缘定位更准噪声严重时反而要降到0.5以下让邻域信息主导分类决策。本文还有配套的精品资源点击获取