CN104881862A - 一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用 - Google Patents

Abstract

本发明属于医学图像处理及应用领域，提供了一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用。本发明首先使用数字化眼底镜获取筛查人群的眼底图像，再利用非下采样的离散小波变换(UDWT)对图像进行增强；然后视网膜灰度图像局部熵的纹理提取，利用模糊C聚类(FCM)的方法对视网膜血管进行分割；最后对分割后的血管进行骨架化，并计算骨架的拓扑层次，并用本发明的迂曲度计算模型对血管骨架进行迂曲度计算。本发明的方法实施简单，精度可靠，便于临床应用。

Description

一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用

技术领域

本发明涉及医学图像处理及应用领域，具体涉及一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用。

背景技术

典型的视网膜图像包括视盘结构，血管和神经由此进入视网膜、从视盘发出的血管分支、视网膜中央的黄斑。正常的视网膜血管走形直或者轻度弯曲，但是在高血流、血管新生和血管阻塞的病理情况下会扩张和迂曲。非光滑的血管走形被称作迂曲，通常是血管形态的首要改变。临床上，迂曲度是反映糖尿病性视网膜病变严重程度的一个重要指标。曲率是形态学分析的一个重要属性，迂曲度通常用血管的曲率变化进行描述，异常的迂曲度是临床许多疾病的重要标志。曲率的估计对于整体眼底血管形状的分析和迂曲的计算具有重要意义。在临床眼科工作中，眼科医生常用定性指标来描述这种改变，如轻度迂曲、中度、重度和极度迂曲。而这种弯曲程度可以计算出数值，对于许多慢性疾病如糖尿病、高血压、心脑血管疾病的疾病的眼底血管弯曲改变定量测量具有重要的意义。

当前，在迂曲度的计算方面，Capowski等人通过测量直的血管上相对增加长度，用一个柔和的曲线表示血管，测量血管积分曲线在一定空间频率的相对长度增加程度，类似的方法也得到了临床应用。但对于测量相对增加长度的方法而言，仅测量曲线的延伸而不涉及形态学和血液动力学结果，而那些通过积分曲线的方法需要在采样的坐标系内通过任意平滑范式来平滑噪声(Capwski JJ,Kylstra JA,Freedman SF.A numeric index based onspatial frequency for the tortuosity of retinal vessels and its application to plusdisease in retinopathy of prematurity[J].Retina,1995,15:490-500.)。Hart等人通过血管曲线的7个积分自动估计迂曲度，但这种方法在某些肉眼看上去就有迂曲改变的图像上没有能分辨出迂曲度的差异(Hart WE,Goldbaum M,Cote B,et al.Automated measurement of retinal vascular tortuosity[J].Int JMed Informatics,1999,53(2-3):239–252.)。Dougherty和Varro等人通过血管中心线的二阶微分来计算迂曲度(Dougherty G,Johnson MJ,WiersMD.Measurement of retinal vascular tortuosity and its application to retinalpathologies[J].Med Biol Eng Comput.2010；48(1):87-95.)。Grisan等人则提出了一种称为曲率相关测量的替代方法，使用曲率的变化来计算(Grisan E,Foracchia M,Ruggeri A.A novel method for the automatic grading of retinalvessel tortuosity[J].IEEE Trans Med Imaging.2008；27(3):310-9.)。这些研究都需要手工分割血管和放置测量点。Sukkaew等人通过两个连续点之间的两个梯度差的和除以采样间隔来计算整个血管的迂曲度系数(Sukkaew L,Uyyanonvara B,Makhanov SS,et al.Automatic tortuosity-based retinopathy ofprematurity screening system.IEICE Trans Inf and Syst,2008,E91-D(12):2868–2874.)。另外，还有两种迂曲度测量矩阵用于对三维结构的血管迂曲度定量分析，并已经被证实在二维和三维血管图像上有效(JohnsonMJ,Dougherty G.Robust measures of three-dimensional vascular tortuositybased on the minimum curvature of approximating polynomial splines fits to thevessel mid-line.Med Eng Phys,2007,29:677-690.Dougherty G,Johnson MJ.Clinical validation of three-dimensional turtuosity metrics based on theminimum curvature of approximating polynomial splines.Med EngPhys,2008,30:190-198.)。Pal等人则提出一种基于链码的迂曲度计算方法，建立k-曲率的改进算法，主要用于提取链码上的分叉点Pi，k用作曲率计算的Pi的数目，一般根据实验以及不同分辨率设置为4,6,8,10,12。K值越小，噪声越多；k值越大又会低估曲率值。传统的迂曲度计算方式τ＝L/C，L为曲线长度，C为曲线端点的直线距离。但这种方法不能判断曲线内部的凹凸程度，因此重新定义迂曲度计算公式为：

$\mathrm{\τ}=\frac{n_{\mathrm{if}}-1}{n_{\mathrm{if}}}\frac{1}{L}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\κ}(P_i,k)$

式中nif表示血管中心线链码中扭转点的数量，L表示弧的长度，这个算法通过对血管上每个像素的曲率求和以及考虑了每个血管亚支的扭转点的数量。除以弧的长度后可以作为迂曲度密度参数，来比较不同长度血管的迂曲度，因此能够避免血管树的分割结果影响(Pal S,Bhowmic P.Estimation of discrete curvature based on chain code pairing and digitalstraightness[C].IEEE Trans ICIP,2009.1097–1100.)。由于正常视网膜血管在视网膜上走行过程中有一定的生理弧度变化，而上述基于曲率的迂曲度表征方法没有能考虑血管的正常弧度，对计算的结果临床意义产生一定的影响。

发明内容

本发明的目的在于提供一种基于眼底镜图像的视网膜血管迂曲度计算方法。

本发明结合视网膜血管在图像上可能存在正常生理弧度这一解剖特点，没有采用传统的曲率计算方法，本发明提出了一种非监督学习的分割方法，并对分割后的视网膜血管骨架进行了拓扑学分层，发明了一种基于随机过程的方法来对迂曲程度进行建模的方法，并对该方法的临床结果进行了验证。

本发明先获取视网膜血管的分割图像，再对视网膜血管进行拓扑分层，利用我们发明的迂曲度计算模型计算拓扑层次上的血管迂曲度。

本发明的技术方案是：经过绿色通道、灰度化等预处理操作后，利用UDWT方法对获得的视网膜图像进行增强，然后视网膜灰度图像局部熵的纹理提取，利用FCM方法对视网膜血管进行分割，对分割后的血管进行骨架化，并计算骨架的拓扑层次，并用我们开发出迂曲度计算模型对血管骨架进行迂曲度计算。

本发明提供了一种基于眼底镜图像的视网膜血管迂曲度计算方法，包括以下步骤：

步骤一：首先使用数字化眼底镜获取筛查人群的眼底图像；

步骤二：利用非下采样的离散小波变换(Undecimated discrete wavelettransform,UDWT)对图像进行增强；

步骤三：然后视网膜灰度图像局部熵的纹理提取；

步骤四：然后利用模糊C聚类(Fuzzy C mean clustering,FCM)的方法对视网膜血管进行分割；

步骤五：对分割后的血管进行骨架化；

步骤六：计算骨架的拓扑层次；

步骤七：利用以下迂曲度计算模型对血管骨架每个层次进行迂曲度计算：

根据血管曲线的每个像素点的X,Y坐标作为输入值，利用协方差矩阵的迹和曲线的像素点个数对血管迂曲度进行建模；

眼底血管待测曲线的每个像素点X,Y坐标为两个随机变量X和Y，根据协方差(记为Cov(X,Y))定义：

Cov(X,Y)＝E{[X-E(X)][Y-E(Y)]}

其中，E(X)＝(X₁+X₂+...Xn)/n；E(Y)＝(Y₁+Y₂+...+Y_n)/n

$\left[\begin{array}{c}{C}_{11}{C}_{21}\\ C_{12}{C}_{22}\end{array}\right]$ 称为r v(X,Y)的协方差矩阵，这里X，Y分别对应于待测血管曲线上的每个像素的横纵坐标值；

迂曲度Tur定义:

Tur＝n²/(C₁₁+C₂₂)

该式中n为待测曲线上的像素点个数，即骨架化后图中像素点个数n，C₁₁+C₂₂为该图骨架化后所有像素点X,Y坐标值的协方差矩阵的迹。

步骤二中的UDWT，可参见文献Starck J L,Fadili J,Murtagh F.Theundecimated wavelet decomposition and its reconstruction[J].IEEE TransImage Process.2007,16(2):297-309.

步骤四中的FCM，可参考文献Yang M S,Hu Y J,Lin K C,et al.Segmentation techniques for tissue differentiation in MRI of ophthalmologyusing fuzzy clustering algorithms[J].Magn Reson Imaging.2002,20(2):173-179.

进一步地，本发明还提供了上述的基于眼底镜图像的视网膜血管迂曲度计算方法在糖尿病早期计算机筛查中的应用。

所述的应用，是指构建糖尿病早期计算机筛查模型、或糖尿病早期计算机筛查系统。

所述的应用包括：步骤二、三、四、五中，利用UDWT和视网膜灰度图像局部熵进行纹理提取，并结合模糊C聚类的算法进行自动视网膜血管的分割方法。

步骤六中，提出的视网膜血管骨架的拓扑层次定义方法。

步骤七中，提出的迂曲度模型对每个层次的血管支进行迂曲度计算，并发现糖尿病早期高层次的血管支发生迂曲。

综上所述，本发明的计算方法具有如下优点：

1.迂曲度是一种无单位的指标，不会受图像的大小，旋转等影响。

2.本发明的方法可以自动化的提取视网膜血管，不需要样本学习。

3.本发明的方法计算简便，运算时间快，敏感度高，在糖尿病患者出现眼底严重病变前发现血管的异常，适合糖尿病筛查。

附图说明

图1为本发明算法的眼底图像的分割结果，其中上排为原始眼底图像，下排为上排对应的分割结果；

图2为本发明提出的视网膜血管树的拓扑层次示意图，其中A用不同形状表示，B用不同数字表示；

□表示起始点或终末点；△表示分支点(每种颜色)；○表示两点间血管层次，相同颜色表示隶属于一个血管树同一层次的血管，第一层1条，第二层2条，第三层4条，第四层4条，第五层2条；

图3为视网膜血管的拓扑分层结果，A为待分层的某个视网膜血管树，B为本发明得到的视网膜血管支层次，C为小波层数设置为1、2、3、4、5、6时的图像；

图4为本迂曲度算法计算的不同模型的迂曲度结果(传统的弧弦比均1.4142)；

图5为本发明的迂曲度算法对同一模型进行不同尺度、旋转变换后的计算结果；

图6为正常和糖尿病眼底血管的分层迂曲度折线图，拓扑层次第9层有统计学意义(P<0.05)。

具体实施方式

下面结合附图对本发明作进一步详细说明，但本发明的实施不仅限于此。

实施例1：

本发明方法通过以下步骤实现：

步骤一，基于FCM的视网膜血管提取

首先对数字眼底镜设备(日本TOPCON)采集到的彩色的RGB眼底图像提取绿色通道，然后进行灰度转换后利用UWT变换，增强后选用第4层的重构图像结构中提取像素点的灰度和局部熵特征，利用FCM算法进行自动血管分割，得到如图1所示血管的分割结果。

本实验中采用首先非下采样离散小波变换(UDWT)对二维图像信号进行分解。小波变换通过平移基本小波(母小波)可获得信号的时间信息，而通过缩放小波的尺度可获得信号的频率特性。对基本小波(母小波)的缩放和平移操作是为了计算小波的系数，这些系数代表小波和局部信号之间的相互关系。本实验中，设定小波母函数:ψ(t)∈L²(R)，其傅里叶变换Ψ(ω)满足如下的允许条件：

$C_{\mathrm{\&Psi;}}=0_R\frac{{\left|\mathrm{\&Psi;}\left(\mathrm{\&omega;}\right)\right|}^2}{\left|\mathrm{\&omega;}\right|}\mathrm{d\&omega;}<\&infin;$

根据视网膜血管纹理自身的特点，本实验中主要用二维离散小波，将图像信号分解成沿行和列的一维问题来处理。离散小波变换针对尺度参数a、平移参数b进行离散化,离散小波变换(DWT)的定义为：

任意L²(R)空间中的x(t)DWT为：其中，获得信号的时间信息，而通过缩放小波的尺度可获得信号的频率特性。对基本小波(母小波)的缩放和平移操作是为了计算小波的系数，这些系数代表小波和局部信号之间的相互关系。二维离散小波可以提取眼底图像的低频系数、水平高频系数、垂直高频系数，对角高频系数，但在传统的离散小波变换在每次计算小波系数的过程中会丢弃前一步的尺度和小波结果，形成类金字塔样分解。而UWT在计算DWT过程中，使用一个N×N环形平移的矩阵TX来代替原来的X矩阵。如果X矩阵是一维的矢量X＝[X₀；X₁；...；XN_-2；XN_-1]，则

TX＝[X_n-1；X₀；X₁；…；X_n-2],T^-1X＝[X₁；X₂；…；X_n-2；X₀]

然后DWT计算一次X，再计算一次TX，并将各自得到的DWT系数进行融合来形成UWT的小波系数【D.B.Percival,A.T.Walden.Wavelet Methods for Time SeriesAnalysis.Cambridge University Press,Cambridge(2000)】。然后，通过对小波系数的重构得到不同尺度下的视网膜图像。

本实验中，通过UWT得到了视网膜图像不同尺度的小波系数，然后根据不同的目的对小波系数进行相应的处理，然后再对处理后的小波系数进行图像重构从而完成图像的处理。重构后的眼底图像就得到了增强。视网膜血管与视网膜图像背景的对比度明显增强，图像清晰(图1)。

本实验选择第四层的UWT重构后图像进行下一步的特征提取。

本方法采用基于模糊C-均值聚类分类器的视网膜血管像素分类。模糊聚类分析就是把一个模糊的不确定的没有类别标记的样本按照某种准则划分为若干个子群集，让相似的样本尽可能归为一类，而把不相似的样本划分到不同的类中。是一种无监督模式识别。模糊聚类算法是一种基于函数最优方法的聚类算法，需要使用微积分计算技术求最优代价函数，在基于概率算法的聚类方法中使用概率密度函数。

模糊集合最根本的方法是其所包含元素及其相应的隶属度用函数表示出来，模糊集合的数学定义如下:

称A是给定论域U上的一个模糊子集的含义是：对任何u∈U都对应一个唯一的实数μA(u)∈[0，1]，有序对集合{u，μA(u)|u∈U}均可作为U上的模糊集，μA(u)表示u对A的隶属度。这表明构造了一个映射：

μ_A(u):U→[0,1],u→μ_A(u)

这个映射称为模糊集A的隶属度函数。一个定义在空间X＝{x}上的隶属度函数就定义了一个模糊集合A，或者叫定义在论域X＝{x}上的模糊子集对于有限个对象x₁，x₂，……，x_n模糊集合可以表示为：

$\underset{~}{A}=\left\{({\mathrm{\&mu;}}_A\left(x_i\right),x_i)\right|x_i\&Element;X\}.$

模糊C-均值聚类(FCM)的核心思想是寻求合适的隶属度和聚类中心，使得聚类内代价函数的方差和迭代误差最小，代价函数的值是像素到聚类中心L²测度的加权累积和，表达形式为:

$J(U,V)=\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ik}}^m{d}_{\mathrm{ik}}^2(x_k,u_i)$

J(U，V)表示区域的像素到聚类中心加权距离的平方和。u_ik表示x_i属于第k类区域的隶属度。m为模糊指数，m∈[1,∞]，表示控制分类矩阵U的模糊程度，是隶属度的加权指数，决定分类结果的模糊程度，其值越大分类越模糊。同时满足约束条件:求解此类极小值的问题，可以通过拉格朗日乘数法来求解。首先，由于矩阵u中各列独立，则有:

$\min \left\{J_m(U,V)\right\}=\min \left\{\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(u_{\mathrm{ik}}\right)}^m{\left(d_{\mathrm{ik}}\right)}^2\right\}=\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\min \left\{\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(u_{\mathrm{ik}}\right)}^m{\left(d_{\mathrm{ik}}\right)}^2\right\}$

求解极小值问题则转化为求解的极小值，且在约束条件 $\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ik}}=1,$ 假设拉格朗日函数F为: $F=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(u_{\mathrm{ik}}\right)}^m{\left(d_{\mathrm{ik}}\right)}^2+\mathrm{\&lambda;}(1-\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ik}}).$ 通过求解偏导数可得出，在约束条件下，其隶属函数和聚类中心的更新依赖于以下两个函数：

$u_{\mathrm{ij}}=\frac{1}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(\frac{d_{\mathrm{ij}}}{d_{\mathrm{kj}}}\right)}^{2/(m-1)}},v_i=\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ik}}^m{x}_k}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ik}}^m}$

由这两个公式得出最佳模糊分类矩阵和聚类中心，这一过程可以通过迭代方法实现^[15]。

本实验过程中，设定FCM的类别数为2，隶属度的加权指数m取2。最大迭代次数为100次，判别误差为1e-5，提取眼底图像中属于血管的像素点。

本发明的方法具有较高的灵敏度(0.82)，算法运行迭代至函数收敛次数平均为50次左右，实现分割的平均时间13.4s。

步骤二，视网膜血管的骨架化和拓扑分层

如图2所示，本发明首先按照图像4邻域法则，识别出眼底血管网络上的特征点如分叉点和终末点。

然后根据建立的特征点，得到两个特征点之间的曲线部分，并根据图2的层次标定方法，对眼底血管骨架网络进行重新标定，标定的结果如图3所示。

本方法提出了一种定义血管树的分支的编号方法，该理论定义所有的子结点都是由根节点派生出来的。结点是按由上往下，由根结点到子节点编号的，首先将根编号发出的分支编号为1；然后根节点向下发出的子节点的分支编号+1，以此类推。(图2)

步骤三，视网膜分层的迂曲度计算

本发明根据血管曲线的每个像素点的X,Y坐标作为输入值，利用协方差矩阵的迹和曲线的像素点个数对血管迂曲度进行建模。

在本发明中，定义眼底血管待测曲线的每个像素点X,Y坐标为两个随机变量X和Y，根据协方差(记为Cov(X,Y))定义：

Cov(X,Y)＝E{[X-E(X)][Y-E(Y)]}

其中，E(X)＝(X₁+X₂+...Xn)/n；E(Y)＝(Y₁+Y₂+...+Y_n)/n

$\left[\begin{array}{c}{C}_{11}{C}_{21}\\ C_{12}{C}_{22}\end{array}\right]$ 称为r v(X,Y)的协方差矩阵，这里X，Y分别对应于待测血管曲线上的每个像素的横纵坐标值。

本发明中定义迂曲度Tur:

Tur＝n²/(C₁₁+C₂₂)

该式中n为待测曲线上的像素点个数，即骨架化后图中像素点个数n，C₁₁+C₂₂为该图骨架化后所有像素点X,Y坐标值的协方差矩阵的迹。

为了验证我们的算法，我们构建出不同的迂曲模型，如图4所示，一峰一谷系列中，峰顶偏右，测得迂曲度为9.4298，峰顶偏左，测得迂曲度为9.3806，峰顶偏高，测得迂曲度为9.4671；两峰一谷系列中，左锋偏左且左右两峰高度差较大，测得迂曲度为9.7881，左锋偏中且左右两峰高度差较大，测得迂曲度为9.7532，左右两峰齐高，测得迂曲度为10.8645，左右峰高度相差较小，测得迂曲度为10.5903；一峰两谷系列中，两谷深度相差较小，测得迂曲度为10.8904，两谷深度相差较大，测得迂曲度为10.4339，两峰两谷系列中，先谷后峰排列的，测得迂曲度为10.3460，先峰后谷排列的，测得迂曲度为10.7003；三峰一谷，测得迂曲度为10.0837；三峰两谷系列中，第三峰较高，测得迂曲度为10.8641，三峰齐高，测得迂曲度为11.2867，第一峰较高，测得迂曲度为11.1562。

如图5所示，取一近半圆曲线的原图，将其分别旋转0°、15°、90°，分别计算其迂曲度，三者均为14.11627；一中部凹陷的曲线，测得迂曲为15.19111，一左部凹陷的曲线，测得迂曲度为14.40097，一右部凹陷的曲线，测得迂曲为14.62633；将一近半圆曲线放大1.2倍，测得迂曲度为14.4027，将一右侧凹陷的曲线放大1.2倍，测得迂曲度为15.50687，将一右侧凹陷的曲线旋转90°，测得迂曲度为15.19111；取一近半圆曲线的原图的截取部分，测得迂曲度为11.1497，将原图两端点距离拉长弧长不变，测得迂曲度为11.7299，将原图两端点距离缩短弧长不变，测得迂曲度为32.8246。

以上结果说明本发明的迂曲度算法能够较好地表征出模型中的迂曲情况。

实施例2：运用本发明方法测量早期糖尿病患者视网膜血管改变

以10幅正常视网膜图像和10幅糖尿病视网膜图像为实验对象，图像有临床眼底镜工作站(南通大学附属医院)提供。图像来自于20例临床确诊为2型糖尿病患者，诊断标准为《1997年美国糖尿病协会提出的糖尿病诊断标准》1、有糖尿病症状，并且任意血糖≧11.1mmol/L。2、空腹血糖≧7.0mmol/L。3、糖耐量试验2小时血糖≧11.1mmol/L。OGTT仍按世界卫生组织的要求进行。由于有10例合并高血压，其眼底照片被排除。最终纳入10幅图像，其中男性6例，女性4例，平均年龄为51.2岁，10例眼底图像按照国际临床DR的严重程度分级和我国糖尿病视网膜病变的分类方法，其中7例被诊断为无明显糖尿病视网膜病变(Diabetic retinopathy,DR)，另外3例被诊断出轻度非增生性DR(有微血管瘤，单纯型I期)。从每一幅图像中选择一段血管枝进行新迂曲度计算。先得到视网膜血管骨架化后的矩阵，然后利用本发明的算法可以实现自动分层标记。在实际血管树图像处理中，会经常发现出血管现交叉情况，对于这类图像，我们设计了一种半自动的血管分层编程方法。首先将眼底血管图像进行分解，其中每个部分内部只有分叉类型，另外图中保留起始端的部分血管枝，便于最后重新组合。再对分解的图进行组合，将其中血管支数据整合成一组数据，再进行使用。

如图6所示，从折线图中可以看出糖尿病眼底血管支的迂曲度起伏不定，最后趋于增大，而正常血管支的迂曲度更加平稳，波动较小。

可见，本发明方法发现早期糖尿病患者在高的视网膜层次上迂曲度提高。

以上已对本发明创造的较佳实施例进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明创造精神的前提下还可作出种种的等同的变型或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。

Claims (2)

1.一种基于眼底镜图像的视网膜血管迂曲度计算方法，包括以下步骤：

步骤一：首先使用数字化眼底镜获取筛查人群的眼底图像；

步骤二：利用非下采样的离散小波变换对图像进行增强；

步骤三：然后视网膜灰度图像局部熵的纹理提取；

步骤四：然后利用模糊C聚类的方法对视网膜血管进行分割；

步骤五：对分割后的血管进行骨架化；

步骤六：计算骨架的拓扑层次；

步骤七：利用以下迂曲度计算模型对血管骨架每个层次进行迂曲度计算：

根据血管曲线的每个像素点的X,Y坐标作为输入值，利用协方差矩阵的迹和曲线的像素点个数对血管迂曲度进行建模；

眼底血管待测曲线的每个像素点X,Y坐标为两个随机变量X和Y，根据协方差记为Cov(X,Y)定义：

Cov(X,Y)＝E{[X-E(X)][Y-E(Y)]}

其中，E(X)＝(X₁+X₂+...Xn)/n；E(Y)＝(Y₁+Y₂+...+Y_n)/n

$\left[\begin{array}{c}{C}_{11}{C}_{21}\\ C_{12}{C}_{22}\end{array}\right]$ 称为rv(X,Y)的协方差矩阵，这里X，Y分别对应于待测血管曲线上的每个像素的横纵坐标值；

迂曲度Tur定义:

Tur＝n²/(C₁₁+C₂₂)

该式中n为待测曲线上的像素点个数，即骨架化后图中像素点个数n，C₁₁+C₂₂为该图骨架化后所有像素点X,Y坐标值的协方差矩阵的迹。

2.一种如权利要求1所述的基于眼底镜图像的视网膜血管迂曲度计算方法在糖尿病早期计算机筛查模型或系统中的应用。