简介这是一份面向医学图像处理学习者的MATLAB视网膜血管分割项目可用于糖尿病视网膜病变、高血压视网膜病变等眼科疾病的辅助分析与早期筛查。项目覆盖完整的血管分割流程图像预处理滤波、对比度增强、光照校正、特征提取Canny边缘检测、Gabor纹理分析、分割算法以及形态学后处理并额外加入图像配准模块帮助对齐不同时间或设备采集的视网膜图像。资源共21个文件以16个m脚本为核心包含预处理、特征匹配、血管连通域处理等功能函数另附tif与png格式的视网膜样本图以及1份PDF文档压缩包仅1.85MB轻量易部署。已有959人学习适合具备一定MATLAB基础、希望深入实践医学图像分割与配准的读者。通过研读源码可以完整掌握从算法原理到代码实现的细节并基于这些模块搭建属于自己的视网膜图像分析流程。1. 视网膜图像血管分割的起点先看清难点在哪拿到一张眼底彩照大多数第一次接触这个题目的人会习惯性先转灰度、再做个大津阈值。我最早做的时候也是这么干的结果分割出来的图像里血管断成一截截能和背景噪声混成一片。视网膜血管分割不是普通的二值化问题它难在血管本身粗细不均、与背景对比度低而且受光照偏斜和视盘边界的影响很大。简单阈值只能把暗斑和粗血管一起捞出来细血管在灰度直方图里几乎和背景重叠根本分不开。这个题目之所以常出现在 matlab 图像处理大作业里是因为它足够考验预处理、滤波和形态学操作的综合运用又不需要特别昂贵的实验环境。本文把这套流程拆成预处理、主干分割、后处理和指标验证四段每一步都给出可以直接在 matlab 命令行里跑的代码你照着搭一遍就能看到从彩照到二值血管图的全过程。2. 绿色通道与形态学顶帽把血管从背景里捞出来2.1 为什么偏偏选绿色通道视网膜彩照是 RGB 三通道直接转灰度是把三个通道的信息平均起来但这个平均会把对比度稀释掉。观察红、绿、蓝三个分量红色通道里血管和背景都偏亮几乎看不出差异蓝色通道噪声最重边缘抖动明显只有绿色通道里血管呈现深色、背景相对明亮两者灰度差最大。这不是玄学眼底血管里的血红蛋白主要吸收绿光所以在绿通道里血管天然更暗。常见的做法是直接取img(:,:,2)后续所有处理都只针对这一个矩阵。这段代码把三通道拆开并显示对比img imread(retina.png); % 原始眼底彩照 R img(:,:,1); G img(:,:,2); B img(:,:,3); montage({R, G, B}); % 并列显示三个通道 title(Red / Green / Blue channel);取绿色通道后下一步不是急着分割而是先做对比度增强。眼底照片经常因为光照不均匀导致一侧偏亮、一侧偏暗直接分割会在暗区产生大片伪影。我一般用adapthisteq做限制对比度自适应直方图均衡化它把图像分成若干小窗口在每个窗口内独立做直方图均衡再通过插值消除块间边界。这样做的好处是暗区的细节能被局部拉伸亮区又不会过曝。% CLAHE 增强关键参数是窗口大小和裁剪限幅 G_enhanced adapthisteq(G, NumTiles, [8 8], ClipLimit, 0.02);这里的NumTiles指把图像分成 8x8 的块块太小噪声放大明显块太大又退化成全局均衡化一般取 [8 8] 或 [16 16] 比较稳。ClipLimit是直方图裁剪阈值默认 0.01 偏保守眼底图像我常用 0.02既能压住背景噪声又能把细血管稍微抬高一点。如果你的图像噪声特别多回落到 0.01 即可如果血管整体偏淡可以往 0.03 方向调整。2.2 顶帽变换滤掉不均匀背景增强之后血管和背景的对比度已经好很多但还有一层障碍背景本身的灰度不是平的视盘周围偏亮、边缘偏暗这种低频分量会给后续阈值分割带来很大麻烦。这里用形态学顶帽变换来去掉背景。顶帽的定义是原图减去开运算结果开运算是先腐蚀后膨胀能后保留比结构元更大的亮区结构把细小的暗血管当噪声滤掉。原图减掉这个结果剩下的就是比结构元小的暗结构正好是血管。se strel(disk, 12); % 圆盘结构元半径要大于血管最大宽度 background imopen(G_enhanced, se); % 估计背景光照 tophat G_enhanced - background; % 顶帽结果血管被保留背景被拉平结构元半径的选择直接影响顶帽效果。半径不能比血管宽否则开运算会把血管也当成背景的一部分减掉。视网膜血管最粗的在视盘附近大概十几个像素我用的是半径 12 的圆盘如果你处理的图像分辨率更高按比例放大到 15 到 20 之间。观察tophat的直方图如果背景的峰值明显集中在 0 附近说明结构元选对了。2.3 对比度增强与顶帽的先后顺序预处理顺序不能乱。先做 CLAHE 再做顶帽和先顶帽再 CLAHE 的结果差别很大。我一般固定为绿色通道提取 → CLAHE → 顶帽变换。原因是顶帽要求背景是缓变的CLAHE 先把光照不均匀压平顶帽才能更干净地分离出血管反过来先顶帽再 CLAHE顶帽阶段会因为原始光照太强而把部分暗血管误判成背景。你可以把两个顺序都跑一遍在后续分割时看血管连续性先 CLAHE 后顶帽的方案细血管断裂明显更少。3. 匹配滤波响应与自适应阈值血管分割的主干算法3.1 匹配滤波为什么适合管状结构血管在局部区域可以看成一段近似直线、剖面呈高斯形的暗色结构。匹配滤波的思路是构造一个与血管剖面形状匹配的二维高斯核把滤波核旋转到不同方向分别和图像做卷积取所有方向响应的最大值。这样无论血管走向是水平、垂直还是斜向总有一个方向的滤波核能对齐它响应值最大而背景中的孤立噪点因为不具备管状轮廓在所有方向上都不会产生一致的强响应。这个方法是 Chaudhuri 提出的经典做法到现在依然是传统血管分割的骨干算法。支持向量机、深度网络等方法当然能在这个任务上做到更高的准确率但匹配滤波不需要训练数据参数少运行速度快在 matlab 里十几行就能实现作为入门和基线非常合适。等你把整套流程跑通了再换深度学习方法时也能知道传统方法的上限在哪、哪些错误类型是它固有的。3.2 构造多方向高斯核高斯核的长度要覆盖血管剖面sigma 控制滤波核的宽度。血管最细的只有 2 到 3 个像素最粗的超过 10 个像素所以 sigma 取 2.5 左右比较折中。核长度取len round(3 * sigma)就够太长了计算量大而且多个方向卷积后边缘伪影重。function kernel vessel_kernel(length_half, sigma, angle_deg) % 生成旋转后的高斯核用于单方向血管匹配 angle angle_deg * pi / 180; [X, Y] meshgrid(-length_half:length_half, -length_half:length_half); Xr X * cos(angle) Y * sin(angle); Yr -X * sin(angle) Y * cos(angle); kernel exp(-Xr.^2 / (2 * sigma^2)); % 只在血管延伸方向附近保留高斯剖面 kernel(abs(Yr) length_half * 0.5) 0; % 零均值化消除平坦区域直流响应 kernel kernel - mean(kernel(:)); end代码里Xr是沿着血管方向的分量Yr是垂直血管方向的分量。高斯分布只在垂直方向上有衰减沿血管方向保持不变这正好模拟理想血管段。abs(Yr) length_half * 0.5这段限制了核的纵向范围避免相邻背景被卷进来。最后必须做零均值化否则图像中平坦区域也会产生非零响应阈值分割时会引入大片假阳性。3.3 多方向响应融合方向步长取 15 度从 0 度到 165 度共 12 个方向即可血管走向不会超出 180 度范围。理论上方向越密匹配越准但 12 个方向已经能达到血管分割的基本要求再多方向只是增加计算时间。tophat double(tophat); % imfilter 对 double 类型处理更稳定 max_response zeros(size(tophat)); for deg 0:15:165 k vessel_kernel(8, 2.5, deg); resp imfilter(tophat, k, same, replicate); max_response max(max_response, resp); end这里用max而不是求平均是因为每个像素只需要响应最大的那个方向如果求平均细血管只在少数方向有强响应会被其它方向的低响应拉低。imfilter的边界选项用replicate它把边缘像素向外复制避免边界处因补零产生暗带。如果你的图像边缘血管比较丰富可以考虑裁剪掉边缘几个像素再做后续处理。3.4 自适应阈值分割匹配滤波响应图里血管区域的值显著高于背景但不同图片的整体灰度范围不一样用固定阈值很脆弱。自适应阈值以响应图的均值和标准差为参照设定阈值为mean k * std。k 越大分割结果越保守只保留最明显的血管细血管丢失多k 越小召回越高但噪声也会混进来。我从经验上先取 k 0.5 观察效果然后按应用场景调整。thresh mean(max_response(:)) 0.5 * std(max_response(:)); mask max_response thresh;如果你希望算法更自动可以用 matlab 的graythresh对响应图做 Otsu 分割。但 Otsu 假设图像是双峰分布匹配滤波响应图不是标准双峰效果常常不如均值加标准差稳健。k 的调节手感比较直观分割结果里血管完整但背景点多增大 k血管断裂很多减小 k。一般来说 k 在 0.2 到 1.0 之间都能找到可用值。4. 后处理与骨架化形态学操作把分割结果修干净4.1 连通域面积滤波去掉孤立噪声阈值分割后的二值图里通常有两类噪声零星分布的孤立亮点以及视盘边缘的不规则块状伪影。前者面积很小只有几个像素后者面积虽然大但形态和血管完全不同。最有效的方法是连通域分析计算每个白色区域的像素个数把小于阈值的区域直接删除。matlab 里bwareaopen就是干这个的。mask_clean bwareaopen(mask, 100); % 删除面积小于100像素的连通域100 像素这个阈值要参考图像分辨率。DRIVE 数据集是 565x584血管主干动辄上千像素小分支也有几十像素100 能保住小分支同时滤掉大部分噪点。如果图像是 2000x2000 级别的眼底相机原图这个阈值要按面积比例放大我一般取 400 到 600 之间。这个操作不会改善血管断裂问题它只负责清掉孤立点。4.2 用形态学闭运算连接断裂处眼底图像里最细的血管只有 1 到 2 个像素宽经过匹配滤波和阈值化后容易断成虚线。闭运算先膨胀后腐蚀能填平小裂缝同时不会像单纯膨胀那样明显改变血管粗细。这里不能用太大的结构元否则会把邻近的细血管黏成一片。se_close strel(disk, 2); mask_connected imclose(mask_clean, se_close);结构元半径 2 只能修复相距最多几个像素的断裂。如果断裂很严重可以先用imdilate把血管整体加粗闭运算后再用imerode缩回原来的粗细本质上是对闭运算的结果做一次形状放缩。但这种操作有风险容易把血管之间的狭小间隙也填上造成血管粘连在血管密集区域尤其明显。4.3 骨取从血管区域到血管骨架血管分割的最终结果在很多应用里不只要一个区域掩码还要血管中心线比如计算动静脉直径比、追踪血管路径。骨架化是提取中心线的标准操作matlab 里bwmorph一行就能完成。skeleton bwmorph(mask_connected, skel, Inf);skel方法反复删除边界像素直到剩下单像素宽的中心线Inf表示一直操作到不能再删为止。骨架结果会有不少毛刺这是血管分叉处和边缘不平整造成的。毛刺会影响后续血管长度测量和分叉点检测可以再做一次端点修剪找出骨架的端点反向追踪删除长度小于 20 像素的短枝。branchpoints bwmorph(skeleton, branchpoints); endpoints bwmorph(skeleton, endpoints); % 从分支点出发找到最近的端点用距离变换修剪短枝 D bwdistgeodesic(skeleton, branchpoints); short_twig D 20 endpoints; skeleton_clean skeleton; skeleton_clean(short_twig) 0;这段代码里bwdistgeodesic计算骨架上的测地距离以分支点为起点向外延伸距离不超过 20 的端点对应的那段路径就是短枝。修剪后骨架保留主干和较长分支后续分叉点检测会更可靠。4.4 后处理顺序的常见误用有一个常见错误是先把骨架提取出来再做连通域过滤这样小噪点会变成小段骨架线在形态上更像血管过滤难度反而更大。我的固定顺序是面积过滤 → 闭运算 → 骨架化 → 短枝修剪。每一步只解决一个问题不要试图用一个大结构元的闭运算同时完成去噪和连接那样会让血管扭曲。处理后的结果肉眼看起来应该是单像素的连续曲线分叉处有少量毛刺可以接受。5. 在 DRIVE 数据上算指标验证分割效果5.1 预测结果和标注怎么对比血管分割的效果不能只看感觉需要量化指标。DRIVE 数据集是视网膜血管分割最常用的基准有 40 张眼底图其中 20 张训练、20 张测试测试集附带人工标注。把上面流程跑出的二值图和标注图对齐逐像素统计四个值真阳性 TP预测为血管且标注为血管、假阳性 FP预测为血管但标注为背景、假阴性 FN预测为背景但标注为血管、真阴性 TN。这些统计量放在混淆矩阵里看。下面的代码假设pred是你的分割结果label是人工标注二者都是逻辑型二值图尺寸已经通过imresize对齐TP sum(pred(:) label(:)); FP sum(pred(:) ~label(:)); FN sum(~pred(:) label(:)); TN sum(~pred(:) ~label(:));这四个值就是后续所有指标的基础。有的资料会把血管标注的粗细规范到统一宽度再对比DRIVE 提供的是手工细标注直接用即可。如果你的标注是从别处生成的注意标注里是否包含视盘区域如果包含报告指标时应说明视盘区域是否被排除。5.2 灵敏度、特异度与准确率的解释从混淆矩阵可以算出一组指标它们从不同角度刻画分割质量指标公式含义准确率 Accuracy(TPTN)/(TPFPFNTN)全部像素中判对的比例灵敏度 SensitivityTP/(TPFN)血管像素中被找出的比例越高漏检越少特异度 SpecificityTN/(TNFP)背景像素中被正确排除的比例越高误报越少Dice 系数2TP/(2TPFPFN)预测与标注的空间重叠程度accuracy (TP TN) / (TP FP FN TN); sensitivity TP / (TP FN); specificity TN / (TN FP); dice 2 * TP / (2 * TP FP FN);灵敏度对细血管的丢失特别敏感如果结果里细血管断得厉害灵敏度会明显往下掉但准确率未必变因为细血管在全体像素里占比太小。所以我通常不只看准确率而是同时盯灵敏度和 Dice。特异性在血管分割里一般很高常常在 0.97 以上它主要反映背景噪声是否被压住。5.3 从分割结果反向调节参数指标不是跑完就完事的要根据指标反推哪一步参数要动。比如灵敏度低、特异度正常说明漏检多应该把自适应阈值的 k 调低或者把匹配滤波的 sigma 调大以捕捉更宽的血管。如果特异度低、灵敏度正常说明噪声被当成了血管需要加大bwareaopen的面积阈值或检查顶帽变换的结构元是否太小。每次只动一个参数记录下来结果变化这样能建立起参数和指标之间的对应关系。5.4 DRIVE 许可与使用提醒DRIVE 的图片来自荷兰糖尿病视网膜病变筛查项目学术研究可以免费使用但要注意它的许可协议禁止在未授权的情况下公开传播原始图片。你在博客或课程报告中展示结果时贴分割结果图比直接贴原图更安全。测试集的标注只用于最终验证不要把它反过来调参与阈值否则指标会虚高真正部署到新图像上时性能会明显回落。如果项目要求更细的血管提取精度那就要从匹配滤波切到 Frangi 滤波加连通域追踪或者直接上 U-Net 这类深度网络但传统方法跑出的这套指标刚好可以作为基线参考。本文还有配套的精品资源点击获取