A Robust Method for Automated Detection and Grading
of Retinal Venous Beading in Diabetes Mellitus
Xu LeiZhang HengyiYu YajunZheng Xiaoxiang
(Dept.of Life Science & Biomedical Engineering, Zhejiang University,Hang Zhou 310027)
AbstractDiabetes mellitus is a high prevalence disease in the world and retinal vessel changes such as venous beading are common findings in patients suffering from long-standing diabetes. Three methods(1.Combined use of the SED, ROSE, Modified LoG algorithm; 2.Gaussian fitting-based method; 3.Gradient-based method) for detection venous caliber were used in the first processing stage, then the magnitude of fast Fourier transform of variation function of venous caliber was used as a detective index of venous beading. The reliability of the results, justified by their low variability, makes feasible the standardization of quantitative parameters to be used both in the diagnosis and in the prognosis of diabetes, and also allows prospective studies based upon them.
Key wordsDiabetes mellitusVenous beadingImage processingGaussian fitting method
1引言
在西方,糖尿病是严重危害人类生命的主要疾病之一,其发病率与死亡率都很高,并且还是失明的一个最常见原因。静脉扩张呈串珠状是糖尿病性视网膜血管病变发展的一个显著特征。因此,该症状严重程度的定量描述,对糖尿病的诊断及预后估计具有十分重要的临床意义。为得到准确性、客观性,及可重复性均能保证的VBD大小、自动、精确的眼底血管管径变化检测算法的提出将是其必要前提。
传统的边沿提取方法分为:梯度法、零交叉边沿检测(LoG算子)法、Canny算子边沿检测法、拟合法、数学形态学方法及动态规划方法等。采用梯度法(Sobel算子等)提取边沿时,对具有理想边沿特性的图像结果是好的(血管边沿为非理想边沿),但同时对噪声也加以放大,且数字图像的非连续性得以加强。而基于2(GI(i,j))过零点的边沿检测算法在选用小参数σ时,虽边沿位置的精度较高,但同时也检测出过多的细小变化,并形成许多不该有的闭合环,易与圆点状出血斑和微血管瘤相混淆。Ulupinar F[1]针对零交叉并不总是对应事实上的边沿,也并不总是准确指示边沿位置的问题(,对边沿进行了修正。Canny J[2]提出了适用于任意形状边沿提取的最优算子(以|(GI(i,j))|局部极大值点作为边沿点)的思想。但由于使用的是高斯滤波器,依然存在着边沿位置不准的问题。利用形态学梯度算子[3]则具有方法简单,速度快(可用专用硬件实现)的优点。但仅适于只含椒盐噪声的图像。Chauduri等[4]采用12个不同方向的高斯模板对眼底图像进行匹配滤波以增强血管部分。但是该方法对血管中轴反光亢进的眼底图像不能很好的加以增强。再加上它使用固定尺寸(模板尺寸:13×9)的高斯模板,所以对管径变化显著及血管高度扭曲的图像效果不佳。若采用不同尺寸的高斯模板将使得计算量聚增。因此,尽管该方法能够在低对比度眼底图像中检测到连续的血管,但只适合于配有专用图像处理硬件的计算机。采用最优滤波策略(直方图方法)或最小费用法(动态规划法)受眼底图像背景不均及眼底生理复杂性因素的影响,处理结果也不好。
根据眼底血管分布本身具有的特点(血管从视神经乳头出发,以放射状扩展至眼底周边部),采用旋转式扫描法进行血管边沿提取及VBD计算,概括为五个步骤:(1)图像采集及前处理;(2)基于SED、ROSE、修正LoG算法融合的血管边沿提取;(3)血管中心线抽取及短枝去除;(4)血管跟踪及管径(直接法、拟合法、梯度法)测量;(5)对管径变化曲线进行FFT及VBD计算。
2方法与步骤
2.1图像采集及前处理
用5%托品酰胺对糖尿病患者的眼睛进行充分散瞳,瞳孔直径达7 mm以上时进行眼底摄影(柯达彩卷)。以视神经乳头及视网膜血管的鼻上、鼻下、颞上、颞下分枝血管的二级分枝处为中心拍摄(视角30°)五张照片,其中视神经乳头图像用作图像拼接的基准。使用激光扫描仪(Uniscan 4C,600 dpi)将眼底照片数字化后(WINDOWS DIB格式)进行平滑滤波(高斯滤波,中值滤波等),消除噪声以提高测量精度。因血管具有良好的形态学特征,采用数学形态学的滚动球算法去除背景信号,然后进行图像增强[5]以提高血管的对比度。
2.2基于SED、ROSE、修正LoG算法融合的血管边沿提取
血管边沿点用公式表示为:LR(i,j)=Rose(Dilation(Bin(SED)(I(i,j))Tem)))·Sgn(SED(I(i,j)))·(LoG(I(i,j)))其中I(i,j)为输入图像。非荧光图像且顺时针扫描时有:LR(i,j)>0代表右边沿点,LR(i,j)<0代表左边沿点。若为荧光图像则反之。Dilation()为膨胀算子,以提高血管粗骨架图像(Bin(SED)(I(i,j))Tem))的连通性(血管对比度差,或者静脉隐匿-Salus, Gunn征均会导致血管间断)。Bin()为二值化算子,Sgn()为符号算子,Tem(Template)表示SED算法[6]中先负后正的模板(宽13Pixels)。表示卷积。ROSE(f)表示以8个不同方向(0°,26.6°,45°,63.4°,90°,116.6°,135°,153.4°)线型结构元(长为18Pixels)所组成的旋转结构元ROSE(Rotating structuring element)作为结构元,对图像f进行开运算[7]。目的在于保留线状的血管部分,去除不需的细节噪声(非血管部分)。
2.3血管中心线抽取及短枝去除
对血管粗骨架图进行十字型中值滤波,以平滑血管边界,避免细化后短枝产生。血管中心线图像用公式表示为:Thin(i,j)=SPTA(MFilter(Dilation(Bin(SED)(I(i,j))Tem)))),MFilter()表示十字型中值滤波,SPTA[8]表示安全点细化算法。采用宽度优先的算法对血管中心线图像进行预跟踪,去除分枝长度小于16点的短枝血管。
2.4血管跟踪及管径测量
根据眼底血管分布的特点,采用视神经乳头区域作为血管跟踪的起始区域(也可人工指定跟踪的起点及结束范围),并在血管分叉处或者动静脉交叉点处跟踪与原跟踪曲线在管径、血管走向、灰度等最为相似的后继血管中心线。每跟踪血管中心线N点(点间间隔为两个像素,取N=16,便于FFT实现)作为一个采样片断。采用三种方法测得延着血管中心线的管径变化曲线。
2.4.1直接测量法采用跟踪血管中心线,同时在当前血管中心线走向的垂直方向(BB′)及垂直方向的附近(±45°,即DD′~EE′内)寻找(血管粗骨架内)连接左边沿点、当前跟踪点、右边沿点这三点距离最短的边线(若不能找到基于LoG算子过零点的边沿点,则以SED算法得到的局部梯度幅度最大值点作为边沿点)。以该连线长度作为当前跟踪点处血管的管径。
2.4.2基于高斯拟合的管径测量方法[4]血管径向灰度模型为:
其中:Ibo=cx+d为背景灰度。
目标函数为:
其中:X=(A,K,μ1,μ2,σ1,σ2,
