简介：基于Hessian矩阵增强的心血管分割是医学图像分析领域的重要课题，这份代码资源面向从事医学影像处理、计算机辅助诊断的研究者与学生，针对血管细长且高对比度结构难以自动提取的痛点，提供一套可运行的MATLAB实现方案。压缩包内共9个文件，以7个m脚本为主，涵盖Hessian矩阵构建、二阶梯度计算、方向估计与Frangi滤波增强等关键步骤；另有1个c文件用于底层图像卷积加速，1个txt文档对代码结构进行说明，整体仅7KB，轻量而便于研读。已有952人学习下载，适合希望掌握基于Hessian矩阵的血管分割原理及代码实现的入门至进阶用户。通过阅读源码，可深入理解Hessian特征值在细长结构增强中的判别作用，并利用阈值筛选、区域生长等后处理方式优化分割效果，为冠状动脉分割研究或课程设计提供直接参考。
1. 血管不是边缘，切不干净才是常态
做过冠状动脉造影分割的人，都会遇到同一个尴尬：血管和骨骼、背景噪声在灰度上高度重叠，你用Canny或Sobel检测出来的边缘断成一截一截，之后无论怎么连都像在拼残图。反直觉的结论是：血管的本质不是“边缘”，而是“管状结构”，想一次提取出完整血管树，需要定位图像局部二阶导数的主方向。基于Hessian矩阵增强的心血管分割，正是利用这一几何性质，将细长血管与斑状噪声、大尺度软组织分离。下面从Hessian特征值分析拆解Frangi滤波器的Matlab实现，从尺度参数到后续分割参数，给出可复现的步骤和调参经验。这套方法适合刚接触医学图像分析的开发者，也适合已有分割基础但想把血管增强做实的人。
2. Hessian矩阵为什么能“看见”血管
2.1 数字图像里的二阶导数怎么算
在连续函数f(x,y)中，Hessian矩阵是二阶偏导组成的对称矩阵，刻画了该点的曲率。数字图像是离散采样，不能直接求导，通常做法是先做高斯平滑，再对高斯核求导得到卷积模板，用卷积完成空间微分。项目里的Hessian2D.m正是这个流程，一个典型的实现如下：
function [Ixx, Ixy, Iyy] = Hessian2D(I, sigma) % sigma : 高斯尺度，与关注的血管直径相关 X = -round(3*sigma):round(3*sigma); G = exp(-X.^2 / (2*sigma^2)); G = G / sum(G); gx = -X / sigma^2 .* G; % 高斯一阶导 gxx = (X.^2 / sigma^4 - 1/sigma^2) .* G; % 高斯二阶导 Ixx = imfilter(imfilter(I, gxx), G', 'replicate'); Ixy = imfilter(imfilter(I, gx), gx', 'replicate'); Iyy = imfilter(imfilter(I, G), gxx', 'replicate'); end这个函数把高斯平滑和偏导计算合并到一块，避免分两步操作带来的额外平滑误差。imfilter的replicate边界填充是为了防止图像边缘因为补零而出现黑色边框。模板半径取3*sigma，是因为高斯核在3倍标准差之外权重已经小到忽略不计。若图像中血管直径只有2~3像素，sigma取1就够；若需要观察粗分支，sigma可能需要到6甚至10。要特别注意，sigma不是越大越好，过大的平滑会把相邻两条血管融成一条。
2.2 特征值分解把局部结构分类
每个像素点得到的2x2 Hessian矩阵是对称的，一定能对角化，分解得到两个特征值 λ1、λ2，并约定|λ1| <= |λ2|。两个特征向量分别指向灰度曲率最小和最大的方向：血管轴向的灰度变化缓慢，因此对应 λ1 接近0；血管法向灰度变化剧烈，因此 λ2 绝对值很大。下表是血管分割中常见的二维局部结构判断规则：
| 特征值条件 | 局部结构 | 造影图中的典型含义 |
|---|---|---|
| λ1≈0，λ2<0 | 亮管状 | 显影血管 |
| λ1≈0，λ2>0 | 暗管状 | 暗血管或管状伪影 |
| λ1<0，λ2<0 | 亮球形/斑块 | 显影剂聚集、钙化点 |
| λ1>0，λ2>0 | 暗球形 | 空腔、低密度阴影 |
| λ1·λ2<0 | 鞍形 | 血管交叉口、不规则背景 |
在心血管造影这类图像中，血管比背景亮，所以最关注“λ1≈0且λ2<0”这一类。eig2image.m就是用来解析计算特征值的，等效公式如下：
s = sqrt(((Ixx - Iyy) / 2).^2 + Ixy.^2); lambda1 = (Ixx + Iyy) / 2 - s; lambda2 = (Ixx + Iyy) / 2 + s;这里s是半轴差，lambda1和lambda2并没有按绝对值大小排序，所以很多实现在计算响应函数前会做一次符号判断。特征值的绝对值大小并不直接说明像素一定是血管，它只负责提供“局部几何形状”的证据，后续 Frangi 响应函数才把这些特征值组合成可用的管状响应。
2.3 梯度响应为什么不适合直接分割血管
梯度即一阶导数，它对任何灰度变化都产生响应。血管边缘会产生两条平行亮线，血管内部反而是平坦区域，因此用梯度幅值阈值得到的是血管“轮廓”而不是“区域”。如果再用形态学填充轮廓内部，又很容易把钙化点一并填进去，后续也难以区分相邻并行血管。Hessian 特征值的优势在于它同时考察两个正交方向的灰度变化率。血管的轴向变化率接近0，法向变化率很大，这种差异是管状结构独有的组合，斑块和噪声很难同时满足。所以在深度学习流行之前，基于Hessian的血管增强一直是医学图像分割的基线方案，即便现在也常被用作预处理的强特征层。
2.4 特征向量方向与血管方向的关系
既然 λ1 对应曲率最小的方向，那这个特征方向就近似是血管的轴向；λ2 对应的特征方向则是跨血管的法向。HessianAng.m做的就是在这个分解基础上计算主方向角度，并把角度转换为图像坐标系下的向量场。这个方向场在后续断点连接中非常有用，后面会在第5章专门说明。需要注意的是，在血管分叉或交叉区域，特征值不再呈现理想的“λ1≈0且λ2<0”，两个特征值的绝对值可能同时变大，方向角也会剧烈跳变。因此不能把Hessian特征方向直接当成可信的中心线方向，而要在增强结果上先过滤管状结构，再用方向场做局部修复。
3. Frangi滤波器的响应函数与多尺度实现
3.1 为什么不用λ2直接做增强
有人会把 λ2 的绝对值直接当作血管响应，结果发现血管边缘响应很高，血管中心反而可能低，而且背景中亮度突变的地方也会出现高响应。这是因为 λ2 只描述了法向曲率的大小，没有排除“边缘”和“斑块”。Frangi 滤波器在1998年提出用组合指标解决这个问题。对二维图像，典型形式可以写成：
V = 0，当 λ1 > 0 或 λ2 > 0 V = exp(-(λ1/λ2)^2 / (2β^2)) * (1 - exp(-(λ1^2 + λ2^2) / (2γ^2)))第一项exp(-(λ1/λ2)^2 / (2β^2))在管状处 λ1/λ2 接近0，该项接近1；在斑块处 λ1≈λ2，该项快速衰减。第二项1 - exp(-(λ1^2+λ2^2) / (2γ^2))是结构抑制项，用来压低平坦背景中的微小波动。β和γ是超参数：β控制“管状判定”的严格程度，γ控制对噪声的容忍度。在FrangiFilter2D.m里默认值通常分别取0.5和15，对绝大多数造影图像能直接使用。
3.2 多尺度是处理粗细血管并存的关键
冠状动脉有主干、分支和末梢，直径可以相差四五倍。固定一个高斯尺度只能增强近似宽度的血管：sigma偏小，主干内部灰度变化平缓，无法形成足够大的二阶响应；sigma偏大，细末梢被高斯平滑抹掉。多尺度做法是使用一组sigma分别计算Hessian和响应，然后对每个像素取所有尺度中的响应最大值。这样粗细血管都在最终增强图中保留下来，同时也会保留一些低置信度的伪结构，需要后续掩膜筛选。
尺度组的选择直接影响分割质量。比较稳妥的经验是：用等差数列，步长0.5，从0.5到6，即scales = 0.5:0.5:6。如果图像分辨率很高、血管主干占40像素以上，可以把上限提高到10。尽量不用等比数列如[1 2 4 8]，因为粗血管区间间隔太大，响应图上容易出现同一根血管不同段的亮度接缝。计算量方面，512×512的图像跑10个尺度在普通PC上约为几百毫秒，可接受。
3.3 可运行的Matlab主循环
结合项目里的Hessian2D.m、eig2image.m，可以写成一个独立函数：
function [enhanced, bestScale] = frangi_2d(I, scales, beta, gamma) % I : 灰度图，double类型 % scales : 高斯尺度向量，例如 0.5:0.5:6 % beta : 管状判定参数，默认0.5 % gamma : 噪声抑制参数，默认15 if nargin < 4, gamma = 15; end if nargin < 3, beta = 0.5; end I = double(I); enhanced = zeros(size(I)); bestScale = zeros(size(I)); for sigma = scales [Ixx, Ixy, Iyy] = Hessian2D(I, sigma); [lambda1, lambda2] = eig2image(Ixx, Ixy, Iyy); % 只保留亮血管：lambda1 < lambda2 < 0 mask = (lambda2 < 0) & (lambda1 < 0); Rb = lambda1 ./ (lambda2 + eps); S = sqrt(lambda1.^2 + lambda2.^2); resp = zeros(size(I)); resp(mask) = exp(-(Rb(mask).^2) / (2*beta^2)) .* ... (1 - exp(-(S(mask).^2) / (2*gamma^2))); updated = resp > enhanced; enhanced(updated) = resp(updated); bestScale(updated) = sigma; end end调用方式：
scales = 0.5:0.5:6; [enh, scaleMap] = frangi_2d(img, scales, 0.5, 15);代码里的eps加在分母上，防止 λ2 为0时除零。mask同时要求两个特征值为负，这是针对亮血管的设定；如果处理的是暗血管，只需对原始图像取负，或者反过来判断正特征值。bestScale保存每个像素的最佳尺度，后续可以用它近似血管直径，这在形态学处理中非常省事。如果直接调用FrangiFilter2D.m这种封装版本，它内部一般也会计算最大响应和尺度索引，但很多实现不会把尺度图返回出来，建议自己改一版保留。
3.4 多尺度参数选择的常见误区
初学常犯的问题是尺度数量太少，比如设[1 4 8]，结果主干增强不连贯。另一个误区是只看增强图数值最大，而没有保存对应的尺度，导致调参时无法定位到底是哪个sigma贡献了响应。还有个容易忽略的细节：当sigma超过15时，卷积模板长度超过90像素，imfilter计算耗时明显上升，但实际新增信息很少，一般没必要设到这么大。最后提醒一句，尺度组上限应参照图像中最大血管的直径，不是越大越好；如果图像里主动脉和冠状动脉同时出现，主动脉壁的曲率会被当成血管增强出来，需要后续用尺度图筛掉过粗结构。
4. 从增强响应到心血管掩膜
4.1 阈值不是随便设的
Frangi增强图的大部分像素值集中在0到0.3，血管中心可到0.9，背景噪声和软组织响应往往低于0.2。直接取0.2作为全局阈值能保留主要分支，但末梢细血管响应很低，容易被切掉。更稳的做法是先用一个 sigma 为0.8的高斯平滑一下增强图，然后按分位数设阈值：
enhSmooth = imgaussian(enh, 0.8); thr = prctile(enhSmooth(:), 95); mask = enhSmooth > thr;分位阈值取95%意味着只保留最亮的5%像素，这是一个起点。如果后续发现主干断裂，把分位数降到90；如果噪声斑块太多，升到98。平滑增强图是为了避免响应图中的孤立尖峰干扰分位数计算。注意imgaussian.m与Hessian2D.m里的高斯平滑是两个用途：前者用于后处理降噪，后者用于计算偏导。
4.2 清除散斑与连接断口
造影剂不一定在整根血管内充盈，末梢常出现断口。形态学闭运算能连接距离很近的断裂：
mask = imclose(mask, strel('disk', 2)); mask = imfill(mask, 'holes');然后按连通域筛选。对每个连通域计算面积、主轴长度和短轴长度，用轴比和面积过滤散斑：
L = bwlabel(mask); props = regionprops(L, 'Area', 'MajorAxisLength', 'MinorAxisLength'); keep = false(size(L)); for k = 1:numel(props) aratio = props(k).MajorAxisLength / (props(k).MinorAxisLength + eps); if props(k).Area >= 30 && aratio >= 1.5 keep(L == k) = true; end end mask = keep;面积阈值30对应约直径5像素的短分支；轴比1.5可以把圆形斑块排除。如果图像分辨率高，需要按面积阈值 = (目标直径/2)^2*pi重新估计。这里的逻辑是：血管是长条形，轴比必然大于1.5；散斑更接近圆形，轴比接近1。闭运算半径2在多数情况下不会把两条并行血管连起来，但若原图分辨率很高且血管间距离很近，需要把半径降到1。
4.3 验证分割是否精准
有金标准标注时，直接用Dice系数和豪斯多夫距离。没有标注时，主要通过三个现象判断：主支是否完整、末梢是否断裂、心腔或主动脉是否被误分割。心腔呈片状，Frangi响应本身不高，但若增强图像平滑过大，心腔边缘的曲线会被误判成粗血管，这时候需要看bestScale的分量：心腔边缘的响应尺度通常远大于冠状动脉直径，可以按直径阈值剔除。
还可以写一个简单的量化脚本辅助调参：
fprintf('连通域数: %d\n', numel(props)); fprintf('平均轴比: %.2f\n', mean(aratio(keep)));当连通域数突然增多但平均轴比下降时，说明阈值放得过宽，混入了碎片。这个数值反馈比肉眼反复对比要稳定得多。
4.4 与区域生长结合提升连通性
阈值分割得到的mask可以作为初选种子区，再在增强图上做区域生长。具体操作是：取增强图响应最高的5%像素作为种子，在原始灰度或增强图上扩展，要求邻域像素与种子像素的灰度差不超过固定值，同时Frangi响应不降到0。这样能补回部分断口，但也容易产生过分割。如果造影剂分布不均，灰度差很难固定；我一般把bestScale作为约束条件，生长方向上的血管直径不能突变，避免连接到旁边的大血管上去。这里的经验是：Frangi增强图只是中间产物，真正成品化时往往会接入区域生长、水平集或骨架连接。
5. 方向场与尺度图用于断点修复和调参
5.1 HessianAng.m 和 grdAng.m 能提供血管方向
HessianAng.m可以根据Hessian特征向量计算局部主角度，这个角度就是血管轴向的方向角。使用方法：
[ang, ~] = HessianAng(lambda1, lambda2, Ixx, Ixy, Iyy);得到的ang分布在0到π之间。拿到方向角后，断口修复就比形态学闭运算精确：先提取分割结果的骨架，对每个断口端点计算方向，两端方向夹角小于30度时，才沿着方向场做插值连接。grdAng.m这类辅助函数通常用来在梯度方向做角度平滑，避免方向场在分叉处突变。
5.2 尺度图可以替代人工半径参数
多尺度循环中保存的bestScale，本质上对应像素处血管直径的代理。把尺度图与mask相乘，能得到每条血管的估计半径，进而动态设置形态学结构元素半径：
radius = max(1, round(bestScale .* mask + 0.5)); maskClosed = false(size(mask)); for r = 1:max(radius(:)) msk = mask & (radius >= r); maskClosed = maskClosed | imdilate(msk, strel('disk', r)); end这样做的好处是：细血管用小半径膨胀，粗血管用大半径闭运算，不容易把并行细血管黏在一起。能想到这一步的项目不多，但一旦用了，后续主干和支干的连接效果会立刻改善。
5.3 快速验证的三招
第一招：增强前先对原始图做分位数归一化I = (I - min(I(:))) / (max(I(:)) - min(I(:)))，让gamma的阈值在不同图像上保持一致。第二招：看bestScale直方图，如果所有血管区的尺度都落在最小值附近，说明scales下限太低或图像本身没有细血管；如果全落在最大值，说明需要提高上限但也要警惕主动脉壁响应。第三招：用montage({uint8(I), imadjust(enh), imdilate(mask, strel('disk',1))}, 'Size', [1 3])把原图、增强图、掩膜放在同一视图里比对，一次迭代就能看出是增强问题还是分割参数问题。
本文还有配套的精品资源，点击获取
