ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

二维Otsu阈值分割原理与Matlab实现:利用邻域均值抗干扰

二维Otsu阈值分割原理与Matlab实现:利用邻域均值抗干扰 简介一个基于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 二维直方图及四个区域定义灰度级数 L256。统计函数 h(i,j) 表示灰度值为 i 且邻域均值为 j 的像素个数其中 j 来自以该像素为中心的窗口平均。除以像素总数得到联合概率 p(i,j)。二维直方图是一个 L×L 矩阵行列分别对应灰度值和邻域均值。阈值向量 (s,t) 将直方图切成四个区域。区域坐标条件语义Af≤sg≤t背景主体Bfsgt前景主体Cf≤sgt低灰度但邻域亮多为边缘Dfsg≤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 同理。如果 ω_A0则这组阈值没有意义直接跳过。类间分离度常用两维迹的形式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, s1, t1); if w0 1e-12, continue; end w1 rectS(S, s2, t2, 256, 256); if w1 1e-12, continue; end mu_i0 rectS(SI, 1, 1, s1, t1) / w0; mu_j0 rectS(SJ, 1, 1, s1, t1) / w0; mu_i1 rectS(SI, s2, t2, 256, 256) / w1; mu_j1 rectS(SJ, s2, t2, 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×3r2 是5×5。kernel 全部元素为 1/(2r1)^2随窗口大小自适应归一化。imfilter 的边界模式用symmetric而不是默认的零填充这样图像最外圈像素的邻域均值不会被人为压暗避免分割结果产生黑边。sub2ind 那行把两个0~255的灰度映射到1~256的线性索引。accumarray 把所有落在同一格子的像素计数一次完成直方图然后除 numel(I) 得到联合概率。如果图像本身已经是逻辑型或uint16im2double 会把它归一化到0~1再量化回0~255这一步顺序不能反过来否则直接 round(img) 对 uint8 会丢失小数。3.3 rectS 函数为什么需要四个边界钳制S 是二维累积和S(i,j) 表示左上角到 (i,j) 的矩形概率总和。要取任意子矩形标准容斥公式需要四个角。真正容易出错的是边界当 s255 时区域B的下界是 s2257已经超过矩阵维度。rectS 在入口判断 i0i1直接返回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 参数速查表参数默认作用调大时的影响r1邻域窗口半径更平滑目标小则丢边界nargin-函数外是否传r不传时使用默认值传0会报错灰度级256直方图宽度减小到64加速峰位移动小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中值滤波再调用本函数。下面表格给出常用的参数组合。场景r预处理后处理文档/票据1无可选形态学开运算光照不均照片2顶帽变换面积连通域过滤强噪声显微图1medfilt2(3)majority过滤遥感地物2CLAHE边缘refinemedfilt2 来自图像处理工具箱处理后噪声点通常消失但图像会轻微模糊。更激进的做法是先用高斯滤波但高斯滤波会模糊边缘导致二维Otsu得到的前景轮廓向内收缩目视上比中值滤波更明显。所以强噪声优先选 medfilt2。4.4 用分离度比值判断可信度主循环里把每次的 sep 存到 sep_grid(s1,t1)循环结束后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在批量图上基本不会翻车。本文还有配套的精品资源点击获取
返回列表