【图像分割】基于分数阶Hessian滤波与自适应主曲率(APC)分析的视网膜血管分割MATLAB实现方法
✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。
🔥 内容介绍
一、背景
（一）视网膜血管分割的重要性
视网膜血管作为人体唯一可以直接观察到的血管系统，其形态和结构的变化与许多全身性疾病（如糖尿病、高血压等）以及眼部疾病（如青光眼、视网膜病变等）密切相关。通过对视网膜图像中血管的准确分割，医生能够获取血管的管径、分支模式、密度等信息，从而辅助疾病的早期诊断、病情监测以及治疗效果评估。例如，在糖尿病视网膜病变的诊断中，血管的异常扩张、新生血管的出现等特征对于判断病情的发展阶段至关重要。因此，视网膜血管分割是医学图像处理领域的研究热点之一。
（二）传统视网膜血管分割方法的局限性
传统的视网膜血管分割方法包括阈值分割、形态学方法、基于机器学习的方法等。阈值分割方法简单快速，但对于光照不均匀、对比度低的视网膜图像，很难选择合适的阈值，容易导致分割不准确，出现血管断裂或误分割的情况。形态学方法通过形态学运算来增强血管结构，但对于复杂的视网膜血管网络，尤其是在血管交叉和弯曲处，可能会丢失部分细节信息。基于机器学习的方法，如支持向量机（SVM）、神经网络等，虽然在一定程度上提高了分割精度，但需要大量的标注数据进行训练，标注过程不仅耗时费力，而且标注结果的主观性可能影响模型的性能。此外，这些方法往往难以有效处理视网膜图像中的噪声和复杂的背景信息。
（三）分数阶 Hessian 滤波与自适应主曲率分析的优势
分数阶 Hessian 滤波能够更好地捕捉图像中的局部结构信息，相比于整数阶 Hessian 滤波，它具有更强的抗噪声能力和对复杂结构的适应性。自适应主曲率（APC）分析可以根据图像局部特征自动调整分析参数，从而更准确地描述视网膜血管的弯曲和分支特性。将分数阶 Hessian 滤波与自适应主曲率分析相结合，有望克服传统方法的局限性，实现对视网膜血管的高精度分割。
三、自适应主曲率 (APC) 分析原理
（一）主曲率的概念
主曲率是描述曲线或曲面弯曲程度的重要几何量。在二维图像中，对于一条曲线（如视网膜血管），其主曲率可以反映曲线在某一点的弯曲程度。主曲率越大，曲线在该点的弯曲程度越大。通过计算图像中各点的主曲率，可以获取血管的弯曲信息，这对于准确分割血管的分支和弯曲部分非常关键。
（二）自适应主曲率分析
传统的主曲率计算方法通常使用固定的参数，对于不同尺度和形状的血管可能不够准确。自适应主曲率（APC）分析则根据图像的局部特征自动调整计算主曲率的参数。具体来说，APC 分析首先对图像进行多尺度分析，例如使用高斯金字塔对图像进行降采样，得到不同尺度下的图像表示。然后，在每个尺度上计算局部区域的主曲率。通过比较不同尺度下的主曲率信息，自适应地选择最适合当前局部区域的主曲率值。这样可以更准确地描述视网膜血管在不同尺度和形状下的弯曲特性。
（三）APC 在视网膜血管分割中的应用
在视网膜血管分割中，APC 分析可以帮助区分血管的主干和分支，以及识别血管的弯曲和交叉区域。通过将 APC 分析得到的主曲率信息与分数阶 Hessian 滤波增强后的图像相结合，可以进一步细化血管的分割结果。例如，在血管分支处，主曲率会发生明显变化，利用这一特性可以准确地定位分支点，从而提高血管分割的精度。
四、基于分数阶 Hessian 滤波与 APC 分析的视网膜血管分割流程
图像预处理：对输入的视网膜图像进行灰度化、降噪等预处理操作，以提高图像质量，减少噪声对后续处理的影响。常用的降噪方法如高斯滤波、中值滤波等。
分数阶 Hessian 滤波：计算图像的分数阶 Hessian 矩阵，并根据上述血管增强函数得到血管增强图像。这一步增强了视网膜图像中的血管结构，使血管与背景的对比度更加明显。
自适应主曲率分析：对增强后的图像进行多尺度分析，计算每个尺度下的主曲率，并自适应地选择最适合局部区域的主曲率值。通过分析主曲率信息，定位血管的分支点和弯曲区域。
分割与后处理：结合分数阶 Hessian 滤波和 APC 分析的结果，使用合适的分割算法（如阈值分割、区域生长等）对图像进行分割，得到初步的血管分割结果。然后，对分割结果进行后处理，如形态学操作（膨胀、腐蚀等），去除孤立的噪声点，连接断裂的血管，进一步优化分割结果。
通过基于分数阶 Hessian 滤波与自适应主曲率分析的方法，可以有效提高视网膜血管分割的精度和可靠性，为医学诊断提供更准确的视网膜血管信息。这种方法充分利用了分数阶 Hessian 滤波对结构的增强能力和自适应主曲率分析对血管弯曲特性的准确描述，克服了传统方法在处理复杂视网膜图像时的局限性。
⛳️ 运行结果
📣 部分代码
%CLAHE
apc_clahe = adapthisteq(apc_mask,'numTiles',[4 4],'nBins',256);
% Thresholding
thresh = isodata(apc_clahe);
apc_thresh = uint8(im2bw(apc_clahe,thresh ));
%If the ISODATA method fails to produce a suitable threshold, a manually
%selected threshold value may be used instead.
% thresh = 0.15;
%apc_thresh = uint8(im2bw(apc_clahe,thresh));
% Filtering out small segments
seg = bwareaopen(apc_thresh, 50);
% ---- local ISODATA function ----
function [level,MAT,MBT]=isodata(Image)
% Reference :T.W. Ridler, S. Calvard, Picture thresholding using an iterative selection method,
% IEEE Trans. System, Man and Cybernetics, SMC-8 (1978) 630-632.
Image = im2uint8(Image(:));
% STEP 1: Compute mean intensity of image from histogram, set T=mean(I)
[counts,N]=imhist(Image);
counts(1)=0;
i=1;
mu=cumsum(counts);
T(i)=(sum(N.*counts))/mu(end);
T(i)=round(T(i));
% STEP 2: compute Mean above T (MAT) and Mean below T (MBT) using T from
% step 1
mu2=cumsum(counts(1:T(i)));
MBT=sum(N(1:T(i)).*counts(1:T(i)))/mu2(end);
mu3=cumsum(counts(T(i):end));
MAT=sum(N(T(i):end).*counts(T(i):end))/mu3(end);
i=i+1;
% new T = (MAT+MBT)/2
T(i)=round((MAT+MBT)/2);
% STEP 3 to n: repeat step 2 if T(i)~=T(i-1)
while abs(T(i)-T(i-1))>=1
mu2=cumsum(counts(1:T(i)));
MBT=sum(N(1:T(i)).*counts(1:T(i)))/mu2(end);
mu3=cumsum(counts(T(i):end));
MAT=sum(N(T(i):end).*counts(T(i):end))/mu3(end);
i=i+1;
T(i)=round((MAT+MBT)/2);
Threshold=T(i);
end
% Normalize the threshold to the range [i, 1].
level = (Threshold - 1) / (N(end) - 1);
MAT = MAT / (N(end) - 1);
MBT = MBT / (N(end) - 1);
end
end