前阵子有个同学拿着课程作业来找我说程序跑出来的边缘图全是断的Sobel阈值调了半天也没救回来。我打开一看问题根本不是阈值而是他连输入图像都忘了先转double。这种场景我见过太多次——图像处理这门课边缘检测和分割永远是新手的第一个坎但这五个算法Roberts、Prewitt、Sobel、Marr-HildrethCanny加Otsu方法本质上是同一条思路串出来的找灰度突变的位置只是各自用了不同的数学工具和工程策略。这篇文章我想按一条完整的脉络把这五个边缘检测算法和一个分割方法从头到尾拆开讲透每一段都配上可以直接跑的Matlab代码并把背后“为什么这么设计”的逻辑也说明白。适合正在做课程作业、准备考试或者刚接触图像处理想系统理解这块内容的人。代码和思路我都会尽量兼顾完整性和可读性你拿去改改就能用。1. 先搞明白边缘检测到底在检测什么1.1 边缘的本质是灰度突变你可以把一张灰度图像理解成一块高低不平的地形——像素灰度值是海拔边缘就是地形上“从一块平坦区域掉到另一块平坦区域”的那条陡坡。陡坡的位置意味着灰度发生了剧烈变化而这种变化在数学上对应着一个东西梯度。一维情况下灰度突变就是函数值跳变导数很大没变化的地方导数为零。二维图像也一样梯度向量的幅度大代表这个像素附近灰度变化剧烈大概率是边缘。所以边缘检测本质上是梯度检测而所有经典边缘检测器的区别都在于“怎么估计梯度”和“怎么判断哪些梯度值得保留”。后面要讲的Roberts、Prewitt、Sobel都是一阶导数的离散近似计算的是梯度本身Marr-Hildreth用的是二阶导检测的是灰度变化率的突变点二阶导过零Canny则是在梯度基础上叠了非极大值抑制和双阈值把“找梯度”升级成“提取有序的边缘线”。理解了这个框架后面每个算法都是往这个骨架里填肉。1.2 为什么不是“一个函数搞定所有事”很多人会问Matlab里edge函数一行就能出结果为什么还要自己写这些算子和原理因为edge是个黑盒子里面每一步都有工程取舍。你不知道它的平滑参数、不知道它怎么处理边界、更不知道为什么某些边缘它选保留某些它选丢弃。遇到效果不合预期时调参就是瞎猜。另外这五个算法各有性格Sobel抗噪不错但边缘粗Marr-Hildreth封闭性好但对噪声极敏感Canny最精细但参数最多Otsu完全不需要人工阈值但只适合双峰明显的直方图。你只有亲手实现一遍才知道这些性格的来源才谈得上选型。1.3 我这篇文章里代码的组织方式我会按照“原理→公式→Matlab实现→效果点评”的顺序来讲每个算法。代码刻意用手写卷积来实现核心逻辑不直接用现成API除了最后对比时会提到edge这样你既能理解原理又能在作业里向老师展示“我真的会”。统一约定所有输入图像先转灰度转成double类型再处理。放码前先给一个读取和转灰度的基础片段%% 读图与预处理 img imread(test_image.png); if size(img, 3) 3 gray rgb2gray(img); else gray img; end f double(gray); % 卷积运算一定用doubleuint8会溢出截断这一步我会在下面每个算法里复用后面不再重复写读取代码。2. Roberts、Prewitt、Sobel三种一阶梯度算子这三个算子都是“小模板滑过图像”的卷积做法区别只在模板多大、权重怎么分配。2.1 Roberts算子2x2模板的直截了当Roberts是1963年提出的老前辈想法非常朴素——用像素对角方向的灰度差来估计梯度。它用两个2x2模板gx [1 0; 0 -1]; % 45度方向差分 gy [0 1; -1 0]; % 135度方向差分先对图像做卷积得到两个方向偏导梯度幅度就是二者平方和再开方G sqrt(gx_img.^2 gy_img.^2)。优点模板最小、计算最快边缘定位精度高。缺点也很明显2x2模板只用4个像素对噪声毫无抑制能力——单像素的椒盐噪声就能同时激活两个模板输出一个强烈的假边缘。另外它只检测对角线方向水平或垂直边缘的响应并不对称导致同一物体在不同角度下边缘强弱不一致。2.2 Prewitt算子3x3模板的均值意义Prewitt把模板扩到3x3解决了Roberts对方向过度敏感的问题。水平模板检测垂直边缘垂直模板检测水平边缘px [-1 0 1; -1 0 1; -1 0 1]; % 检测垂直边缘 py [-1 -1 -1; 0 0 0; 1 1 1]; % 检测水平边缘拆开看就能发现Prewitt模板其实做了两件事先对图像做3x1的均值平滑同一列取三个值累加再做水平差分。均值平滑能压掉一部分随机噪声所以同样面对带噪图像时Prewitt输出的假边缘会比Roberts少一截。代价是什么均值平滑本身会把边缘“抹平”一点导致检出的边缘比真实边缘略宽——这就是为什么Prewitt边缘图看起来总比Roberts“粗一圈”。边缘定位精度有所牺牲换取噪声抑制能力这是一个典型的工程权衡。2.3 Sobel算子的加权诀窍Sobel算子是Prewitt的升级版模板从全1变成了中心权重2sx [-1 0 1; -2 0 2; -1 0 1]; % 检测垂直边缘 sy [-1 -2 -1; 0 0 0; 1 2 1]; % 检测水平边缘多加的权重2意味着什么中心像素的灰度对梯度贡献更大。这等价于先用一个带权重的平滑核[1;2;1]在垂直方向平滑再做水平方向的中心差分。中心加权的平滑核更接近高斯形状对噪声的抑制效果比Prewitt的均匀均值更好同时保持了边缘定位的锐利度。这也是为什么Sobel在工程上比Prewitt更常用——性价比极高模板只有3x3但效果逼近更复杂的算子。2.4 三种算子的Matlab手写实现%% Roberts roberts_gx conv2(f, gx, same); roberts_gy conv2(f, gy, same); roberts_mag sqrt(roberts_gx.^2 roberts_gy.^2); %% Prewitt prewitt_gx conv2(f, px, same); prewitt_gy conv2(f, py, same); prewitt_mag sqrt(prewitt_gx.^2 prewitt_gy.^2); %% Sobel sobel_gx conv2(f, sx, same); sobel_gy conv2(f, sy, same); sobel_mag sqrt(sobel_gx.^2 sobel_gy.^2); %% 显示直接显示梯度幅度从数据范围归一化到0-255 figure; subplot(2,2,1); imshow(f, []); title(原图); subplot(2,2,2); imshow(roberts_mag, []); title(Roberts梯度幅度); subplot(2,2,3); imshow(prewitt_mag, []); title(Prewitt梯度幅度); subplot(2,2,4); imshow(sobel_mag, []); title(Sobel梯度幅度);注意两点conv2默认在边界补零这在图像边缘区域会产生明显的渐变假边缘——后面我在“常见问题”章节里会单独说边界处理。另外imshow加[]参数是让Matlab自动把最小值映射到0、最大值映射到255避免梯度值范围不一导致图太黑或太白。2.5 用阈值把梯度图变成边缘图上面的梯度幅度图是连续灰度图还不是“边缘图”。要得到二值边缘需要做一次阈值化梯度大于一定阈值的像素标记为边缘。阈值怎么定最简单的是设成最大梯度的某个百分比thresh 0.2 * max(sobel_mag(:)); edge_binary sobel_mag thresh;这个0.2就是“灵敏度旋钮”。调低了边缘多但是噪声也多调高了边缘干净但是会断。没有绝对标准得看你的图和需求。如果你不想手调也可以直接用Matlab内置自动阈值thresh_auto graythresh(sobel_mag / max(sobel_mag(:))); edge_binary sobel_mag thresh_auto * max(sobel_mag(:));graythresh是按Otsu方法计算阈值这个我们第五章节细讲。3. Marr-Hildreth算子用二阶导的零交叉找边缘3.1 为什么还要绕到二阶导一阶导找的是梯度极大的位置但梯度变化是渐变的尤其是边缘带有过渡带或模糊时一阶导的峰值很平缓难以精确定位。二阶导呢峰值点对应的是二阶导的过零点——从正变负或从负变正的那一个像素位置定位比一阶导的峰更精确。Marr-Hildreth也叫LoG算法Laplacian of Gaussian的思路先用高斯核平滑图像消除噪声再对平滑结果求拉普拉斯二阶导。因为线性滤波可以合并实际实现时直接用一个“高斯拉普拉斯”核一次卷积完成。3.2 LoG核的数学形式与生成连续二维LoG函数是LoG(x, y) -1/(πσ^4) * [1 - (x^2y^2)/(2σ^2)] * e^(-(x^2y^2)/(2σ^2))这个算子在中心是正或负看符号约定的碗状外围是一圈反向的环形象地说像个墨西哥草帽。用Matlab生成LoG核function h log_kernel(sigma, n) % n为核尺寸一般取6*sigma取奇 if nargin 2 n 2 * ceil(3 * sigma) 1; end half (n - 1) / 2; [x, y] meshgrid(-half:half, -half:half); r2 x.^2 y.^2; h -(1 / (pi * sigma^4)) * (1 - r2 / (2 * sigma^2)) .* exp(-r2 / (2 * sigma^2)); endsigma决定核的尺寸和抗噪能力。sigma越大平滑越强检出的边缘越“宏观”但定位误差也会变大sigma太小噪声压不住会冒出大量假边缘。3.3 零交叉检测的Matlab实现卷积后用LoG核处理图像得到的是二阶导响应图。边缘在哪在响应从正到负跨越的那个位置。判断方法遍历每个像素检查它4邻域或8邻域内有没有正负号变化H log_kernel(2, 15); LoG_resp conv2(f, H, same); % 零交叉检测 [rows, cols] size(LoG_resp); edge_log zeros(rows, cols); for i 2:rows-1 for j 2:cols-1 patch LoG_resp(i-1:i1, j-1:j1); if max(patch(:)) 0 min(patch(:)) 0 edge_log(i, j) 1; end end end这段循环有O(n²)复杂度图片大时会慢但理解起来最直观。工程优化版本可以用形态学膨胀配合符号比较来做向量化这里不展开了。实话说LoG有一个显著特点输出的边缘是闭合的——它不像Sobel那样只在梯度大的地方出线它会把所有灰度变化的边界都圈出来包括缓慢变化的阴影边界。这在某些场景下是优点在另一些场景下就是灾难因为缓变区域也会被算作边缘图会显得很“脏”。实际用的时候通常还会加一个梯度幅度阈值过滤掉弱响应。4. Canny边缘检测工程界最常用的标准方案Canny在1986年提出三个目标低错误率真实边缘不漏检、非边缘不多检、定位精准检出的边缘像素尽量贴近真实边缘、单像素响应同一条边缘只输出一条线而不是一坨。为了同时满足这三个目标Canny设计了五步流水线。我认为理解Canny的窍门是把它“拆成流水线”每一步解决一个具体问题。4.1 五步流水线第一步高斯平滑去噪。这一步直接影响后续所有结果和LoG一样sigma太小压不住噪声sigma太大会抹掉细边缘。第二步计算梯度幅度和方向。这一步通常用Sobel算子实现得到每个像素的梯度幅度和方向4方向量化水平、垂直、45度、135度。第三步非极大值抑制——把梯度峰值两边不是最大的像素置零让边缘变细。第四步双阈值——设置高阈值TH和低阈值TL把梯度大于TH的像素标记为强边缘介于TL和TH之间的标记为弱边缘小于TL的丢弃。第五步滞后连接——强边缘一定是最终边缘弱边缘只有与强边缘连通才保留孤立弱边缘丢弃。这一步能有效去除“断裂的噪声响应”。4.2 非极大值抑制的细节非极大值抑制是Canny里最容易被忽视的步骤。直觉上讲边缘附近的梯度响应是一个“小山丘”峰值在边缘正中心两侧梯度逐渐降低。如果我们把所有高于某阈值的像素都标成边缘那边缘会是一条宽带。非极大值抑制就是在梯度方向上只保留局部最大的那个像素。实现要点判断一个像素是不是局部最大要在它梯度的方向上看它的邻居。梯度方向垂直于边缘所以要让3x3邻域落在梯度方向上的两个邻居非格点方向都做比较。通常把梯度方向量化为四种角度范围0度、45度、90度、135度然后比较对应位置的梯度幅度。例如梯度方向接近0度水平方向就比较左右两个水平邻居的梯度值。例如量化后方向为0度当前像素值若大于左边邻居且大于右边邻居就保留否则置零。下面给出一个完整但较容易理解的双循环版本注意用向量化优化% 假设 Gmag 是梯度幅度Gdir 是梯度方向弧度 Gdir_q round(Gdir / (pi/4)); % 量化 % 按方向决定比较邻居篇幅有限完整非极大值抑制代码会有点长库函数edge(I,canny)内部已经封装好了。教学实现时可以用简化版只沿梯度方向比较两邻域值的最大值来更新梯度图。4.3 双阈值和滞后连接设高阈值TH低阈值TL通常TL 0.5 × TH。强边缘直接保留弱边缘只有“贴”着强边缘才保留。怎么判断“贴着”常用做法是8邻域递归/广度优先传播从强边缘像素出发向8邻域扩散遇到弱边缘像素就把它提升为强边缘再继续扩散直到没有新像素可提升为止。这个传播过程会让原本断裂的边缘重新连起来——弱边缘如果孤立存在大概率是噪声则没有强邻居支撑删除。4.4 Matlab实现手写版 vs 内置版如果你要手写Canny我建议分文件实现逐步调试主框架如下%% Canny手写主流程 sigma 1.5; G fspecial(gaussian, 2*ceil(3*sigma)1, sigma); f_smooth conv2(f, G, same); [Sx, Sy] gradient(f_smooth); % 数值梯度替代Sobel Gmag sqrt(Sx.^2 Sy.^2); Gdir atan2(Sy, Sx); % 非极大值抑制见上小节 % 双阈值滞后连接 TH 0.3 * max(Gmag(:)); TL 0.1 * max(Gmag(:));如果你就是为了出结果、做对比直接用内置函数就好edge_canny edge(f, canny, [TL TH], sigma);注意这里的[TL TH]如果直接用0到1的范围值Matlab会当成归一化阈值也可以不传就让Matlab自动估计阈值。实际上手写版的价值在于让你能逐行打断点看中间结果真正理解“为什么Canny边缘是连续单像素的”——这一点理解了后面调任何边缘检测参数都有底气。用Canny出的边缘细、连续、定位准这也就是为什么它在工业界地位这么高。5. Otsu方法自动找阈值的图像分割5.1 分割与阈值的连接边缘检测告诉你“哪里是边界”图像分割要告诉你“哪些像素属于同一个区域”。最朴素的区域分割就是灰度阈值分割选一个灰度值T灰度大于T的像素标为前景小于T的标为背景。问题是T怎么选人工调的话不同图片不同光照得反复试Otsu就是为了解决这个“自动选T”的问题。5.2 最大类间方差的思想Otsu的准则非常清晰如果阈值T把灰度分成前景C0和背景C1两类我们希望这两类之间的差异尽可能大、每一类内部尽可能均匀。类间方差的公式为σ²B(T) ω0(T) * ω1(T) * [μ0(T) - μ1(T)]²其中ω0和ω1分别是前背景像素比例μ0和μ1是两类的平均灰度。遍历所有可能的灰度值T0到255让σ²B最大的那个T就是最优阈值。这个准则非常直观如果阈值选在直方图的两个峰中间的“谷底”那么两类均值差距最大、类内方差小类间方差最大如果阈值偏到某个峰的内部两类均值差距变小类间方差自然也变小。所以找最大类间方差等于自动找谷底。5.3 Matlab手写Otsufunction T otsu_threshold(f) % 输入f为double灰度图 counts imhist(uint8(f)); % 256个灰度级的像素计数 total numel(f); p counts / total; % 灰度概率分布 omega0 cumsum(p); mu_cum cumsum((0:255) .* p); mu_total mu_cum(end); % 遍历所有阈值1~255 sigmaB zeros(256, 1); for t 1:255 omega0_t omega0(t); omega1_t 1 - omega0_t; if omega0_t 0 || omega1_t 0 sigmaB(t) 0; continue; end mu0 mu_cum(t) / omega0_t; mu1 (mu_total - mu_cum(t)) / omega1_t; sigmaB(t) omega0_t * omega1_t * (mu0 - mu1)^2; end [~, T] max(sigmaB); % 找最大类间方差对应的灰度 end调用时T otsu_threshold(f); segmented f T;注意imhist(uint8(f))——如果你输入的是doubleMatlab的imhist默认把0~1里的值当成灰度级别直接调用得到的是错误直方图。转成uint8是正确姿势。另外如果图像直方图是单峰比如大部分像素是背景前景极小Otsu的类间方差曲线不会有明显的峰自动选出的阈值常常偏离预期这时要考虑用固定阈值或先做直方图均衡化。5.4 内置函数对比Matlab内置对应的是T_norm graythresh(gray); % 返回0~1的归一化阈值 T_real round(T_norm * 255); bw imbinarize(gray); % 直接得到二值图graythresh统计过程和我上面的代码等效所以你自己写一遍的意义更多是理解而不是超越。不过手写版本会让你真正明白“阈值的取值依据是直方图形状而不是某个固定灰度值”这件事。6. 常见问题与调参技巧实录6.1 噪声是边缘检测的头号敌人Roberts和Marr-Hildreth对噪声极敏感。前者因为模板太小没平滑后者因为拉普拉斯本身就放大噪声。所以实际用的流程永远是先平滑再检测边缘。Canny内部已经做了高斯平滑不需要你再先平滑一次但如果你觉得Canny检出的边缘噪声多优先降低sigma而不是中途再模糊一次——重复模糊会让边缘偏移。6.2 阈值选择的标准做法对于Sobel/Prewitt这类梯度图我通常按梯度最大值的百分比选阈值初始选0.15到0.3之间观察边缘连续性再微调。如果场景复杂、边缘多而杂可以改用自适应阈值比如局部均值加常数替代全局阈值。对于Canny双阈值里高阈值控制边缘的“主干信度”低阈值控制“延伸灵敏度”经验上低阈值设成高阈值的一半左右再慢慢微调。调试时建议实时显示二值图和梯度幅度图对照比只看最终结果更容易判断问题出在阈值还是出在预处理。6.3 Matlab数据类型踩坑清单这部分是实操中遇到最多的坑我直接列出来uint8卷积后结果仍然是uint8负数会截断成0正数超过255截断成255所以运算前一定转double。conv2默认边界补零图像边缘一圈会产生很大伪梯度。用imfilter(f, kernel, replicate)或者symmetric来延伸边界能显著减轻边缘效果。例如prewitt_gx imfilter(f, px, replicate)。graythresh输入图像是uint8或double都可以但返回阈值永远是0到1归一化值。你要实际灰度阈值时别忘了乘以255上面提过这个坑这里再强调一次。显示用imshow(img, [])自动拉伸灰度范围避免因为你转成double后范围不在0到255而显示成一片黑。imhist(uint8(f))和histogram(f(:), 256)在统计灰度分布时结果一致但细节有区别调试直方图时留意。6.4 各算法怎么选型我用一个表来归纳具体场景下的选型建议需求场景推荐算法关键参数注意事项实时性要求高、图像干净Roberts无噪声环境下极易误检教学演示、基础边缘Prewitt无边缘较宽需配合阈值通用场景、默认首选Sobel阈值均衡性好工程常用需要闭合边缘轮廓Marr-Hildrethsigma噪声敏感需先降噪边缘细腻、定位精度要求高Cannysigma、双阈值参数多但效果稳定自动分割而不想调阈值Otsu无直方图双峰时效果好目标与背景灰度重叠很大以上都难—考虑纹理、颜色或深度信息从经验上看如果不知道选哪个默认走SobelCanny的组合——Sobel看整体趋势Canny看细节轮廓两者互补足够应付大多数入门项目。6.5 手写边缘检测的调试心得手写实现时我强烈建议“逐步可视化”。也就是每算完一个中间结果就figure; imshow(...)看一眼。比如写Canny先看平滑后的图再看梯度幅度图再看NMS后的图最后看双阈值连接结果。哪一步不对一眼就知道。我曾调试过一个非极大值抑制的实现发现边缘全错位就是因为量化方向时把atan2结果量化错了象限图一对比就暴露了。这种逐层调试的经验比背十遍公式都有用。7. 写在最后的一点经验回顾整个算法家族你会发现它们本质上都在回答同一个问题如何用离散的像素数据精确地估计连续灰度场中的突变位置。Roberts是最朴素的第一代Prewitt和Sobel是工程化的第二梯队Marr-Hildreth绕到二阶导获得更锐利的定位Canny用非极大值抑制和双阈值把边缘从“点集”变成“线”。而Otsu虽然名义上是分割但它解决的问题和边缘检测殊途同归——都是把图像里的“结构”从“背景”里剥离出来。对我个人来说学这块最大的收获是不要直接信任API也不要恐惧手写。你亲手写一遍Roberts的2x2模板会真正理解为什么它抗噪差你亲手实现一遍Canny的滞后连接会彻底搞懂双阈值的意义。这些理解会在你后续做分割、目标检测、轮廓提取时反复用上。如果你正在写课程报告建议把每种算法的中间结果图都贴出来配上你的阈值设定和调试说明这比只贴一张最终效果图有说服力得多。最后送你一个小技巧跑代码之前先看一眼图像的直方图心里对灰度分布有个概念再选择边缘检测还是Otsu分割——有时候问题根本不是算法而是你一开始就选错了工具。