ARTICLE DETAIL

资讯详情

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

MATLAB手写Canny边缘检测:从高斯滤波到双阈值原理与实现

MATLAB手写Canny边缘检测:从高斯滤波到双阈值原理与实现 简介面向图像处理初学者与MATLAB使用者的Canny边缘检测算子实现资源核心价值在于引导读者以手动编程方式复现经典Canny算法而不仅依赖内置edge函数。Canny算子由John Canny于1986年提出以其低错误率、精准边缘定位和单像素响应成为边缘检测领域的事实标准广泛应用于遥感图像分析、医学影像处理、工业缺陷检测等任务。资源包为zip格式包含2个文件docx算法说明文档与m格式MATLAB源码压缩包整体仅22KB轻巧且无冗余文件。docx文档系统拆解了算法全流程——高斯滤波平滑去噪、Sobel算子计算梯度幅值与方向、非极大值抑制细化边缘、双阈值检测连接强弱边缘并结合实际参数选择给出讲解m脚本按同样顺序逐行实现用户可运行调试亦可修改核尺寸或阈值观察效果。当前已有868人学习对正在准备课程设计、毕业设计或希望打下扎实图像处理基础的读者是一份简洁但信息量充足的参考资料。1. Canny 算子凭什么还能打MATLAB 手写实现的价值做图像处理时我见过太多人一行edge(img,canny)出完图然后换一张带噪声的实拍图就直接翻车。阈值不对、边缘断裂、伪边缘比真边缘还亮。这个 1986 年的算法在今天的工业视觉、医学影像和遥感图里依然是默认基线——因为它的四步结构里每一步都在跟噪声和定位精度做权衡。如果你用的正好是 MATLAB除了内置功能也可以自己实现一遍 Canny才能真正理解高斯滤波、梯度计算、非极大值抑制和双阈值的作用。这份资源包里的canny.m和 docx 说明文档正好带你走一遍这条路。我最初是拿它做 matlab 图像处理大作业后来在一条视觉检测产线上改缺陷定位也靠手写版本把误报率压下来。所以这篇东西对两类人有用一类是学生想知道edge()背后到底做了什么另一类是工程师需要在特定场景下替换梯度算子、固定阈值或插入自定义后处理。2. 四步实现 Canny滤波、梯度、非极大值抑制、滞后连接2.1 高斯滤波平滑噪声但不能牺牲边缘定位Canny 的第一步不是算梯度而是先平滑。原因很简单图像噪声在差分后会变成极大的梯度值如果不先压制噪声后面的非极大值抑制会把噪声当边缘保留下来。高斯滤波用二维高斯核做加权平均离中心越近权重越大这样既平滑又不会像均值滤波那样把边缘步长破坏得太厉害。高斯滤波的唯一参数是标准差sigma核尺寸必须跟着它走。常见做法是取hsize 2 * ceil(3 * sigma) 1原因是高斯分布在 ±3σ 以外的权重已经接近 0。工程上我更倾向于用这种方式而不是imgaussfilt(img, sigma)因为你能明确看到核形状方便在导出报告或做嵌入式移植时替换成定点数。下面这段是标准的滤波入口% 读取图像并转为 double避免后续卷积和差分出现截断 img imread(demo.png); if size(img, 3) 3 img rgb2gray(img); end img double(img); % uint8 下负梯度会溢出必须转换 sigma 1.0; hsize 2 * ceil(3 * sigma) 1; % sigma1.0 时核为 7x7 gkernel fspecial(gaussian, hsize, sigma); img_sm imfilter(img, gkernel, replicate);imfilter第三个参数replicate表示边界像素复制扩展这比默认的补零行为更接近真实场景。补零会在图像四周制造一条暗边随后在梯度阶段被误检成边缘。如果你用imgaussfilt它内部默认处理了边界但可调的只有 sigma无法直接控制核尺寸。需要做定制滤波时还是上面这组代码更稳。2.2 梯度幅值和方向选 Sobel 还是中央差分平滑之后需要计算每个像素的梯度幅值和方向。Sobel 算子包含一个中心差分和一个局部平滑比单纯使用gradient()的中央差分对噪声更鲁棒。在 Canny 的标准实现里Sobel 是默认选项换成 Prewitt 或 Scharr 也可以但差异主要体现在转角处的响应强度。[Gx, Gy] imgradientxy(img_sm, sobel); [Gmag, Gdir] imgradient(Gx, Gy); % Gdir 单位是度范围 [-180, 180]imyradientxy返回的是水平、垂直方向的一阶导数imgradient再合成幅值和角度。角度方向指向梯度变化最快的方向即垂直于边缘的方向后面非极大值抑制就是沿着这个方向找局部最大值。你也可以自己用conv2写 Sobelsobel_x [-1 0 1; -2 0 2; -1 0 1]; sobel_y sobel_x; gx conv2(img_sm, sobel_x, same); gy conv2(img_sm, sobel_y, same);对比一下两种梯度算子算子类型实现方式在 Canny 中的表现Sobel3x3 分离核带一点平滑默认选择抗噪好Prewitt3x3 均匀核更简单边缘稍粗伪响应少Central diffgradient()函数噪声放大明显一般先配合imgaussfilt使用在 MATLAB 里我一般直接用imgradientxy因为conv2要自己处理边界而imgradientxy默认的边界消隐方式更干净。手写版本如果想导出 C 代码再拆成conv2不迟。2.3 非极大值抑制边缘细化算法的实现细节梯度幅值图里一条边缘往往是“亮带”而不是“亮线”因为边缘两侧都会产生梯度响应。非极大值抑制NMS的作用是保留沿梯度方向的局部最大幅值把亮带压成单像素宽。实现时通常把梯度方向量化到四个区域水平-22.5°~22.5°、45°22.5°~67.5°、垂直67.5°~112.5°、135°112.5°~157.5°。然后沿梯度方向比较当前像素和相邻两个像素的幅值。这里给出一个可直接放进canny.m的子函数function nms my_nms(mag, dir) % 非极大值抑制沿梯度方向比较相邻像素幅值 [rows, cols] size(mag); nms zeros(rows, cols); dir mod(dir, 180); % 角度统一切换到 [0, 180) for r 2:rows-1 for c 2:cols-1 d dir(r, c); if d 0 d d 180; end if (d -22.5 d 22.5) || (d 157.5) % 水平梯度方向比较左右邻点 if mag(r, c) mag(r, c-1) mag(r, c) mag(r, c1) nms(r, c) mag(r, c); end elseif d 22.5 d 67.5 % 45 度方向比较右上、左下 if mag(r, c) mag(r-1, c1) mag(r, c) mag(r1, c-1) nms(r, c) mag(r, c); end elseif d 67.5 d 112.5 % 垂直方向比较上下 if mag(r, c) mag(r-1, c) mag(r, c) mag(r1, c) nms(r, c) mag(r, c); end else % 135 度方向比较左上、右下 if mag(r, c) mag(r-1, c-1) mag(r, c) mag(r1, c1) nms(r, c) mag(r, c); end end end end end这个函数只处理内部像素图像最外一圈直接置零。多数应用里边缘不会贴着图像边界出现这个损失可以接受。如果确实需要完整尺寸输出可以先用padarray把原图扩展一圈做完 NMS 再裁掉。2.4 双阈值检测和滞后连接用形态学重建一步到位NMS 之后仍有大量小响应需要用双阈值筛选。高于高阈值的点一定是边缘低于低阈值的点一定不是介于两者之间的点只有与强边缘连通才保留。这就是“滞后连接”。一个避免写循环的常见做法是用形态学重建imreconstruct(strong, weak, 8)强边缘作为种子弱边缘作为掩膜重建后保留下来的弱边缘就是与强边缘连通的点。代码很短% 阈值取梯度最大幅值的比例比例需要根据图像内容调整 highThreshold 0.2 * max(Gmag(:)); lowThreshold 0.4 * highThreshold; strong nms highThreshold; weak nms lowThreshold; % weak 包含 strong作为重建掩膜 % imreconstruct 从 strong 出发在 weak 范围内做 8 邻域连通重建 edges imreconstruct(strong, weak, 8);这里weak必须包含strong因为重建算法要求种子在掩膜内部。网上不少手写实现把双阈值写成edge_map (Gmag low) | (Gmag high)这个写法有几个问题一是括号没有把和|隔开实际逻辑会被优先级影响二是丢失了 NMS 结果导致边缘变粗。正确做法是先对 NMS 结果做阈值再用imreconstruct做连接。如果你关心每一步的中间形态可以把 strong、weak 单独保存后面调试阈值时很有用。3. 封装成可复用的 my_canny.m并与内置 edge() 对比3.1 把四个步骤收进一个函数写项目代码时我习惯把 Canny 封装成函数参数暴露sigma、低阈值比例、高阈值比例方便在批处理脚本里循环调参。下面是一个精简版封装内部调用了上一节的my_nmsfunction [edges, Gmag, Gdir] my_canny(img, sigma, lowRatio, highRatio) % 简易 Canny 实现返回二值边缘和中间变量 if nargin 2, sigma 1.0; end if nargin 3, lowRatio 0.4; end if nargin 4, highRatio 0.2; end if size(img, 3) 3 img rgb2gray(img); end img double(img); % 高斯滤波 hsize 2 * ceil(3 * sigma) 1; gkernel fspecial(gaussian, hsize, sigma); img_sm imfilter(img, gkernel, replicate); % 梯度计算 [Gx, Gy] imgradientxy(img_sm, sobel); [Gmag, Gdir] imgradient(Gx, Gy); % 非极大值抑制 nms my_nms(Gmag, Gdir); % 双阈值与滞后连接 highTh highRatio * max(Gmag(:)); lowTh lowRatio * highTh; strong nms highTh; weak nms lowTh; edges imreconstruct(strong, weak, 8); end调用方式很直观img imread(cameraman.tif); [edges, Gmag, Gdir] my_canny(img, 1.0, 0.4, 0.2); imshow(edges);注意lowRatio和highRatio是比例关系具体数值取决于图像内容。在my_canny里我把中间变量Gmag、Gdir也返回了这样调试时不需要重新计算梯度直接拿来做直方图和可视化。3.2 手写实现和内置 edge()差异在哪MATLAB 内置的edge(img, canny)同样是四步结构但有几个关键差别对比项内置 edge()手写 my_canny梯度算子默认 Sobel内部可调固定调用 imgradientxy非极大值抑制使用线性插值比较四方向量化比较高阈值由梯度累计直方图自动估计手动指定比例滞后连接内部 bwmorph / 邻域逻辑imreconstruct 重建用一句代码就能把两个结果叠在一起看差异edgeBuiltin edge(img, canny); imshowpair(edges, edgeBuiltin, blend);从实际效果看内置版本在弱边缘连接上更连续因为它的 NMS 用插值得到亚像素级比较手写版在斜线上会有轻微断点。但手写版的优势是你能改任意一环。比如产线上印刷图案有方向性我直接把手写版的 Sobel 换成 Scharr边缘的连续性立刻好了很多内置edge()并不提供这个替换入口。3.3 中间结果可视化定位故障在哪一步调 Canny 最怕“边缘不对但不知道哪一步错”。把中间结果用montage排在一起是最直接的诊断手段montage({ mat2gray(img), ... mat2gray(img_sm), ... mat2gray(Gmag), ... mat2gray(nms), ... edges ... });正常情况下Gmag应该看不清明显的方向性纹理nms是细线edges里边缘连续、没有碎片。如果Gmag出现了大量噪点检查sigma是否太小如果nms还是宽条检查角度量化时mod(dir, 180)是否正确如果edges断裂严重把lowRatio调大到 0.5 左右再跑一遍。4. 参数怎么选高斯核、阈值比例与 MATLAB 里的常见坑4.1 sigma 和核尺寸的对应关系高斯滤波的sigma直接影响边缘的定位精度。sigma 越大越能压噪声但强边缘会被磨圆sigma 太小噪声产生的高频梯度会让 NMS 输出一堆碎点。实际处理时可以参考下面这个范围σ核尺寸示例适用场景0.53x3 或 5x5高质量图片纹理细密1.05x5 或 7x7常规照片、文档扫描1.57x7 或 11x11中等噪声如手机室内拍摄2.011x11 或 13x13强噪声低照度监控场景核尺寸必须是奇数这样锚点才能在正中心。如果发现整幅图的边缘都向同一方向偏移多半是核尺寸取了偶数或者imfilter的边界选项没设好而不是 Canny 本身的问题。4.2 双阈值比例用梯度直方图自动估计大多数教程直接给low0.4*high和high0.2*max(Gmag)但碰到边缘密集的图像按最大值取阈值会把真正的弱边缘全部滤掉。更稳的做法是看梯度幅值的累计分布用分位点代替固定比例% Gmag 来自 my_canny 的输出edgesHist 是直方图的 bin 边界 [counts, edgesHist] histcounts(Gmag(:), 100); cum cumsum(counts) / sum(counts); % 取累计分布 80% 处作为高阈值这样只有 20% 的像素可能成为边缘 idx find(cum 0.8, 1, first); highTh edgesHist(idx 1); lowTh 0.4 * highTh;用分位点而不是最大值能避免某个极端高亮像素把高阈值顶上去。比如图像里有一个很强的镜面反光max(Gmag(:))可能比正常边缘高 10 倍按固定比例算出来的高阈值会大得离谱导致真实边缘全部丢失。分位点在批处理不同图片时更鲁棒。4.3 三个容易改错的细节第一uint8运算溢出。imfilter或gradient处理uint8图像时差分结果可能变成负数MATLAB 会截断到 0导致梯度方向直接错误。所以我在 2.1 节一开始就把图像转换成double。这一点对 Canny 是致命的但对普通平滑滤波影响不大所以很多人发现不了。第二角度单位不统一。atan2(Gy, Gx)返回弧度必须先转成度再做mod(dir, 180)。如果直接用弧度做象限判断-3和3在弧度下相差很大量化区间全部错位NMS 后边缘会变成双线。第三和|的优先级。MATLAB 中优先于|所以a b | c会被解析成(a b) | c。双阈值逻辑建议全部加括号例如weak (nms lowTh) (nms highTh)避免后面复习代码时还要猜意图。4.4 分步输出把中间变量写成图片调试 Canny 时我习惯在my_canny.m里临时加几行imwrite把中间变量导出成图片再放到imtool里对比imwrite(mat2gray(img_sm), debug_1_sm.png); imwrite(mat2gray(Gmag), debug_2_grad.png); imwrite(mat2gray(nms), debug_3_nms.png); imwrite(edges, debug_4_edge.png);导出而不是直接用imshow好处是可以切换不同的可视化工具放大看像素值。nms阶段最值得看如果nms图里有断断续续的暗点说明阈值偏低NMS 没有完全抑制噪声如果强边缘周围出现白色光晕说明高斯核太大了边缘定位偏移。5. 把 Canny 的结果用起来边缘方向直方图5.1 从二值边缘提取方向分布Canny 输出的二值边缘除了拿来显示还能继续统计方向信息。在工业视觉里判断一块工件是否放正不需要跑完整的目标检测直接统计边缘方向直方图就能看出主方向。做法是把Gdir和edges对应起来dirEdge Gdir(edges); dirEdge(dirEdge 0) dirEdge(dirEdge 0) 180; % 统一到 0~180 histogram(dirEdge, 36, Normalization, probability); xlabel(梯度方向度); ylabel(概率);梯度方向与边缘本身的方向相差 90°。比如竖直边缘的梯度方向是 0° 或 180°水平边缘的梯度方向是 90°。实际使用时如果看到直方图在某个角度出现明显的峰就说明图像里有大量该方向的边缘结构。5.2 过滤小连通域避免噪声拉平直方图图像上的孤立噪点在被 Canny 检测为边缘后会产生随机方向这些方向会把直方图的主峰拉平。统计前先用bwareaopen去掉小连通域edgesClean bwareaopen(edges, 10); % 去掉面积小于10像素的连通域 dirClean Gdir(edgesClean); histogram(dirClean, 36, Normalization, probability);bwareaopen的第二个参数是面积阈值。10 像素对大多数视觉检测来说足够小不会误伤真实边缘。如果边缘本身很破碎可以把这个值降到 5但代价是噪声方向也会进来。实践中我会先看edges里有多少个面积小于 10 的孤立点数量占比超过 5% 才调整阈值。5.3 计算主方向并做角度校正有了方向直方图就能估算主边缘方向。下面的代码用histcounts取峰值 bin 的中心角度[hCounts, binEdges] histcounts(dirClean, 0:5:180); [~, peakIdx] max(hCounts); dominantAngle binEdges(peakIdx) 2.5; % bin 中心 if dominantAngle 90 dominantAngle dominantAngle - 180; % 转成 -90~90 end fprintf(主边缘方向: %.1f 度\n, dominantAngle);得到dominantAngle后可以直接作为imrotate的反向角度来校正图像。注意这里算出来的是梯度方向边缘方向需要再加 90° 并取模 180否则旋转方向会相反。这个方法在 PCB 定位、木材纹理方向检测这些场景里比跑深度学习要快得多而且每张图只需一次 Canny 再加一次直方图统计。本文还有配套的精品资源点击获取
返回列表