一种糖尿病性眼底静脉串珠状变形程度定量分析方法

http://www.100md.com 《生物医学工程学杂志》 1999年第1期

     作者：许 雷 张恒义 虞亚军 郑筱祥

    单位：许 雷 张恒义 虞亚军 郑筱祥(浙江大学 生命科学与医学工程学系，杭州 310027)

    关键词：糖尿病；串珠状静脉；图像处理；高斯拟合

    生物医学工程学杂志990119 内容摘要 糖尿病性眼底静脉串珠状变形程度VBD(Venous Beading Degree)的定量分析，对糖尿病的无损诊断及预后估计具有十分重要的临床价值。采用三种算法对静脉管径进行自动获取。1.基于SED、ROSE、修正LoG算法^[1]融合的血管边沿提取及管径计算。2.基于高斯拟合的管径求取。3.基于梯度的管径求取。对测得的管径变化曲线经线性插值后进行FFT，由Fourier频谱的大小得到定量描述的VBD。结果表明：运用本文方法所得到的VBD在正常人与糖尿患者之间具有显著性差异，并具有精度高、可重复性好及可自动测量的优点。

    A Robust Method for Automated Detection and Grading

    of Retinal Venous Beading in Diabetes Mellitus

    Xu Lei Zhang Hengyi Yu Yajun Zheng Xiaoxiang

    (Dept.of Life Science & Biomedical Engineering, Zhejiang University,Hang Zhou 310027)

    Abstract Diabetes 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 words Diabetes mellitus Venous beading Image processing Gaussian fitting method

    1 引言

    在西方，糖尿病是严重危害人类生命的主要疾病之一，其发病率与死亡率都很高，并且还是失明的一个最常见原因。静脉扩张呈串珠状是糖尿病性视网膜血管病变发展的一个显著特征。因此，该症状严重程度的定量描述，对糖尿病的诊断及预后估计具有十分重要的临床意义。为得到准确性、客观性，及可重复性均能保证的VBD大小、自动、精确的眼底血管管径变化检测算法的提出将是其必要前提。

    传统的边沿提取方法分为：梯度法、零交叉边沿检测(LoG算子)法、Canny算子边沿检测法、拟合法、数学形态学方法及动态规划方法等。采用梯度法(Sobel算子等)提取边沿时，对具有理想边沿特性的图像结果是好的(血管边沿为非理想边沿)，但同时对噪声也加以放大，且数字图像的非连续性得以加强。而基于²(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] 血管径向灰度模型为：

    其中：I_bo=cx+d为背景灰度。

    目标函数为：

    其中：X＝(A,K,μ₁,μ₂,σ₁,σ₂,c,d)为模型参数；M为径向灰度采集点个数(M=19)；为惩罚因子。

    图1 血管管径的直接测量

    Fig 1 Directly measuring vessel's caliber

    运用优化方法求出拟合参数σ₁，则血管管径近似为：VC≈｜6σ₁｜。此法的优点在于抗噪声能力强。

    2.4.3 基于梯度的管径测量方法 设血管截面近似为圆，则血管边沿点的粗确位置由下式O(i)的局部极大值所在的位置确定^[9]：

    式中：

    其中：I_i为采集到的血管径向灰度值；I^′_i为对I_i求一阶梯度；I^″_i为对I_i求二阶梯度。

    2.5 对管径变化曲线进行FFT及VBD大小计算

    受图像离散化的影响，沿着血管中心线的管径测量间隔及血管径向灰度采样间隔是不等的(1或)。因此需采用线性插值的方法对管径曲线进行重采样，使得采样间隔等距化，然后再进行FFT。设原始管径变化曲线为d′(l);l(l:1～N pixels)是沿着血管中心线的离散采样距离，插值后的管径曲线为：d(n)。其离散Fourier变换为：，则VBD定义为：。显然，VBD越大则静脉串珠状变形程度越大。

    3 结果与讨论

    图2为三幅图像的叠加图像，即：(1)原始眼底图像(Size：512×512)，(2)采用基于SED、ROSE、修正LoG算法融合方法得到的血管边沿图像，(3)血管中心线及血管径向灰度采集线图像。图3为一采集到的血管径向灰度曲线及其基于高斯的血管径向灰度拟合曲线(σ₁=2.665；σ₂=0.597；μ₁=-0.281；μ₂=0.035；c=0.255；d=162.9；A＝17.93；K=32.78)。运用直接法对正常人及糖尿病患者的VBD进行测量后得：正常人NVBD=9.0±3.02，糖尿病患者为DVBD＝30.3±16.83，具有显著性差异，临床意义显著。本文介绍了三种方法，对管径进行了测量，虽测得的管径值并不完全一样，但测得的管径变化却是一致的。且在管径获取的精度上(SNR较高时)，梯度法>直接法>拟合法。在抗噪声能力上，拟合法>直接法>梯度法。所需的处理时间：拟合法>梯度法>直接法。因此，在实际应用中应根据图像的质量、机器的速度等具体条件加以选择。利用本文方法也可对眼底血管的血管迂曲率，动静脉管径比等血形态参数进行自动定量分析，因此本文方法具有很高的临床使用价值。

    图2 血管边沿、血管中心线及血管径向灰度采集线示意图

    Fig 2 Illustration of vessel edges, vessel center lines and sampling lines of vessel intensity profile

    图3 典型的视网膜血管径向灰度曲线及其相应的高斯拟合曲线

    Fig 3 A typical intensity profile of a retinal vessel and its related Gaussian-based fitting curve

    *浙江省自然科学基金资助项目(JD9503)

    参考文献

    [1] Ulupinar F, Gerard Medioni. Refining edges detected by LoG operator. Comput Vis Graph and Image Process, 1990; 51∶275

    [2] John Canny. A computational approach to edge detection. IEEE Trans PAMI-8,1986∶679

    [3] Peng J, Rusch Ph, Herrmann T, et al. Morphological filters and edge detection application to medical imaging. Annual International Conference of the IEEE Engineering in Medicine and Biology Society. 1991; 13(1)∶0251

    [4] Chauduri S, Chatterjee S, Katz N et al. Detection of blood vessels in retinal images using two-dimensiolal matched filters. IEEE Trans Med Imaging, 1989;8∶263

    [5] 许 雷，郑筱祥.一种基于视觉非线性及数学形态学的视网膜血管图像自适应增强方法.中国医疗器械杂志，1997； 21(3)∶144

    [6] Tascini G, Passerini G, Puliti P et al. Retina vascular network recognition. Proc. SPIE, 1993; 1898∶322

    [7] Brett DT, Alan CN. Semi-automatic segmentation of vascular network images using a rotating structuring element (ROSE) with mathematical morphology and dual feature thresholding. IEEE Trans Med Imaging, 1993；12∶385

    [8] Naccache NJ, Shinghal RS. A proposed algorithm for thinning binary patterns. IEEE Trans System, Man, and Cybernetics, 1984,14(3)∶409

    [9] Steven RF, Maryl RJ, Christopher J et al. Automated analysis of coronary arterial morphology in cineangiograms:Geometric and physiologic validation in humans. IEEE Trans Med Imaging, 1989,8∶387

    (收稿：1997-06-13)

    百拇医药网 http://www.100md.com/html/analecta/2003/08/31/73/688.htm