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

一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用
【技术领域】
[0001] 本发明设及医学图像处理及应用领域，具体设及一种基于眼底镜图像的视网膜血 管迂曲度计算方法及其应用。
【背景技术】
[0002] 典型的视网膜图像包括视盘结构，血管和神经由此进入视网膜、从视盘发出的血 管分支、视网膜中央的黄斑。正常的视网膜血管走形直或者轻度弯曲，但是在高血流、血管 新生和血管阻塞的病理情况下会扩张和迂曲。非光滑的血管走形被称作迂曲，通常是血管 形态的首要改变。临床上，迂曲度是反映糖尿病性视网膜病变严重程度的一个重要指标。曲 率是形态学分析的一个重要属性，迂曲度通常用血管的曲率变化进行描述，异常的迂曲度 是临床许多疾病的重要标志。曲率的估计对于整体眼底血管形状的分析和迂曲的计算具有 重要意义。在临床眼科工作中，眼科医生常用定性指标来描述该种改变，如轻度迂曲、中度、 重度和极度迂曲。而该种弯曲程度可W计算出数值，对于许多慢性疾病如糖尿病、高血压、 屯、脑血管疾病的疾病的眼底血管弯曲改变定量测量具有重要的意义。
[0003] 当前，在迂曲度的计算方面，Capowski等人通过测量直的血管上相对增加长 度，用一个柔和的曲线表示血管，测量血管积分曲线在一定空间频率的相对长度增加程 度，类似的方法也得到了临床应用。但对于测量相对增加长度的方法而言，仅测量曲 线的延伸而不设及形态学和血液动力学结果，而那些通过积分曲线的方法需要在采样 的坐标系内通过任意平滑范式来平滑噪声（CapwskiJJ,KylstraJA,化ee血anSF.A numericindexbasedonspatialfrequencyforthetortuosityofretinal vesselsanditsapplicationtoplusdiseaseinretinopathyofprematurity[J]. Retina, 1995, 15:490-500.)。Hart等人通过血管曲线的7个积分自动估计迂曲度， 但该种方法在某些肉眼看上去就有迂曲改变的图像上没有能分辨出迂曲度的差异 (HartWE,GoldbaumM,CoteB,etal.Automatedmeasurementofretinalvascular to;rtuosity[J].IntJMedInformatics, 1999, 53 (2-3): 239 - 252.)。Dou曲erty和 Varro等人通过血管中屯、线的二阶微分来计算迂曲度值ou曲ertyG,^hnsonMJ,Wiers MD.Measurementofretinalvasculartortuosityanditsapplicationtoretinal pathologies[J].MedBiolEngComput. 2010 ;48(1): 87-95.)。Grisan等人则提出 了一种 称为曲率相关测量的替代方法，使用曲率的变化来计算（GrisanE，ForacchiaM，Ruggeri A.Anovelmethodfortheautomaticgradingofretinalvesseltortuosity[J]. IE邸hansMedImaging.2008;27(3):310-9.)。这些研究都需要手工分割血管和放 置测量点。Sukkaew等人通过两个连续点么间的两个梯度差的和除臥采样间隔来计算 整个血管的迂曲度系数（SukkaewL，UyyanonvaraB，MaldianovSS，etal.Automatic tortuosity-basedretinopathyofprematurityscreeningsystem.lEICETransInf andSyst，2008，E91-D(12):2868 - 2874.)。另外，还有两种迂曲度测量矩阵用于对S 维结构的血管迂曲度定量分析，并已经被证实在二维和S维血管图像上有效（Johnson MJ,DoughertyG.民obustmeasuresofthree-dimensionalvasculartortuositybased ontheminimumcurvatureofapproximatingpolynomialsplinesfitstothe vesselmid-line.MedEngPhys,2007,29:677-690.Dou組ertyG,JohnsonMJ.Clinical validationofthree-dimensionalturtuositymetricsbasedontheminimum curvatureofapproximatingpolynomialsplines.MedEngPhys,2008,30:190-198.)。 Pal等人则提出一种基于链码的迂曲度计算方法，建立k-曲率的改进算法，主要用于提取 链码上的分叉点pi，k用作曲率计算的Pi的数目，一般根据实验W及不同分辨率设置为 4,6,8, 10, 12。K值越小，噪声越多；k值越大又会低估曲率值。传统的迂曲度计算方式T =L/C，L为曲线长度，C为曲线端点的直线距离。但该种方法不能判断曲线内部的凹凸程 度，因此重新定义迂曲度计算公式为：
[0004]
[0005] 式中nif表示血管中屯、线链码中扭转点的数量，L表示弧的长度，该个算法通过 对血管上每个像素的曲率求和W及考虑了每个血管亚支的扭转点的数量。除W弧的长度 后可W作为迂曲度密度参数，来比较不同长度血管的迂曲度，因此能够避免血管树的分割 结果影口向（PalS,BhowmicP.Estimationofdiscretecurvaturebasedonchaincode pairinganddigitalstrai曲tness[C].I邸ETransICIP, 2009. 1097 - 1100.)。由于正 常视网膜血管在视网膜上走行过程中有一定的生理弧度变化，而上述基于曲率的迂曲度表 征方法没有能考虑血管的正常弧度，对计算的结果临床意义产生一定的影响。

【发明内容】

[0006] 本发明的目的在于提供一种基于眼底镜图像的视网膜血管迂曲度计算方法。
[0007] 本发明结合视网膜血管在图像上可能存在正常生理弧度该一解剖特点，没有采 用传统的曲率计算方法，本发明提出了一种非监督学习的分割方法，并对分割后的视网膜 血管骨架进行了拓扑学分层，发明了一种基于随机过程的方法来对迂曲程度进行建模的方 法，并对该方法的临床结果进行了验证。
[0008] 本发明先获取视网膜血管的分割图像，再对视网膜血管进行拓扑分层，利用我们 发明的迂曲度计算模型计算拓扑层次上的血管迂曲度。
[0009] 本发明的技术方案是；经过绿色通道、灰度化等预处理操作后，利用UDWT方法对 获得的视网膜图像进行增强，然后视网膜灰度图像局部滴的纹理提取，利用FCM方法对视 网膜血管进行分割，对分割后的血管进行骨架化，并计算骨架的拓扑层次，并用我们开发出 迂曲度计算模型对血管骨架进行迂曲度计算。
[0010] 本发明提供了一种基于眼底镜图像的视网膜血管迂曲度计算方法，包括W下步 骤：
[0011] 步骤一：首先使用数字化眼底镜获取筛查人群的眼底图像；
[001引步骤二：利用非下采样的离散小波变换（Undecimateddiscretewavelet transform,UDWT)对图像进行增强；
[0013]步骤S;然后视网膜灰度图像局部滴的纹理提取；
[0014] 步骤四；然后利用模糊C聚类（FuzzyCmeanclustering,FCM)的方法对视网膜 血管进行分割；
[0015] 步骤五：对分割后的血管进行骨架化；
[0016] 步骤六：计算骨架的拓扑层次；
[0017] 步骤走：利用W下迂曲度计算模型对血管骨架每个层次进行迂曲度计算：
[0018] 根据血管曲线的每个像素点的X，Y坐标作为输入值，利用协方差矩阵的迹和曲线 的像素点个数对血管迂曲度进行建模；
[0019] 眼底血管待测曲线的每个像素点X，Y坐标为两个随机变量X和Y，根据协方差（记 为Cov狂，Y))定义；
[0020] Cov化Y) =E{技-E狂)][Y-E(Y) ]}
[00引]其中，E佩二化+X2+. . .Xn)/n;E(Y)=化+Y2+. . .+YJ/n
[0022]
称为rv〇(，Y)的协方差矩阵，该里X，Y分别对应于待测血管曲线上的每 个像素的横纵坐标值；
[002引 迂曲 度Tur定义：
[0024] Tur = n]/ (〔。+〔22)
[00巧]该式中n为待测曲线上的像素点个数，即骨架化后图中像素点个数n，C11+C22为该 图骨架化后所有像素点X，Y坐标值的协方差矩阵的迹。
[0026]步骤二中的UDWT，可参见文献StarckJIMuda曲F.Ilie undecimatedwaveletdecompositionanditsreconstruction[J].IEEETransImage Process. 2007, 16(2) : 297-309.
[0027]步骤四中的FCM,可参考文献YangMS,HuYJ,LinKC,etal.Se卵entation techniquesfortissuedifferentiationinMRIofophthalmologyusingfuzzy clusteringalgorithms[J].MagnResonImaging. 2002, 20(2):173-179.
[0028] 进一步地，本发明还提供了上述的基于眼底镜图像的视网膜血管迂曲度计算方法 在糖尿病早期计算机筛查中的应用。
[0029] 所述的应用，是指构建糖尿病早期计算机筛查模型、或糖尿病早期计算机筛查系 统。
[0030] 所述的应用包括：步骤二、S、四、五中，利用UDWT和视网膜灰度图像局部滴进行 纹理提取，并结合模糊C聚类的算法进行自动视网膜血管的分割方法。
[0031] 步骤六中，提出的视网膜血管骨架的拓扑层次定义方法。
[0032] 步骤走中，提出的迂曲度模型对每个层次的血管支进行迂曲度计算，并发现糖尿 病早期高层次的血管支发生迂曲。
[0033] 综上所述，本发明的计算方法具有如下优点：
[0034]1.迂曲度是一种无单位的指标，不会受图像的大小，旋转等影响。
[00巧]2.本发明的方法可W自动化的提取视网膜血管，不需要样本学习。
[0036] 3.本发明的方法计算简便，运算时间快，敏感度高，在糖尿病患者出现眼底严重病 变前发现血管的异常，适合糖尿病筛查。
【附图说明】
[0037]图1为本发明算法的眼底图像的分割结果，其中上排为原始眼底图像，下排为上 排对应的分割结果；
[003引图2为本发明提出的视网膜血管树的拓扑层次示意图，其中A用不同形状表示，B用不同数字表示；
[0039] □表示起始点或终末点；A表示分支点（每种颜色）；0表示两点间血管层次，相 同颜色表示隶属于一个血管树同一层次的血管，第一层1条，第二层2条，第=层4条，第四 层4条，第五层2条；
[0040] 图3为视网膜血管的拓扑分层结果，A为待分层的某个视网膜血管树，B为本发明 得到的视网膜血管支层次，C为小波层数设置为1、2、3、4、5、6时的图像；
[0041] 图4为本迂曲度算法计算的不同模型的迂曲度结果（传统的弧弦比均1.4142);
[0042] 图5为本发明的迂曲度算法对同一模型进行不同尺度、旋转变换后的计算结果；
[0043] 图6为正常和糖尿病眼底血管的分层迂曲度折线图，拓扑层次第9层有统计学意 义（P<0. 05)。
【具体实施方式】
[0044] 下面结合附图对本发明作进一步详细说明，但本发明的实施不仅限于此。
[0045] 实施例1 ;
[0046] 本发明方法通过W下步骤实现：
[0047] 步骤一，基于FCM的视网膜血管提取
[004引首先对数字眼底镜设备（日本T0PC0N)采集到的彩色的RGB眼底图像提取绿色通 道，然后进行灰度转换后利用UWT变换，增强后选用第4层的重构图像结构中提取像素点的 灰度和局部滴特征，利用FCM算法进行自动血管分割，得到如图1所示血管的分割结果。
[0049] 本实验中采用首先非下采样离散小波变换扣DWT)对二维图像信号进行分解。小 波变换通过平移基本小波（母小波）可获得信号的时间信息，而通过缩放小波的尺度可获 得信号的频率特性。对基本小波（母小波）的缩放和平移操作是为了计算小波的系数，该 些系数代表小波和局部信号之间的相互关系。本实验中，设定小波母函数：(t)G。(时， 其傅里叶变换W(?)满足如下的允许条件；
[0050]
[0051] 根据视网膜血管纹理自身的特点，本实验中主要用二维离散小波，将图像信号分 解成沿行和列的一维问题来处理。离散小波变换针对尺度参数a、平移参数b进行离散化， 离散小波变换值WT)的定义为：
[005引任意L2佩空间中的X(t)DWT为
获得信号的时间信息，而通过缩放小波的尺度可获得信号的频率特性。对基本小波（母小 波）的缩放和平移操作是为了计算小波的系数，该些系数代表小波和局部信号之间的相互 关系。二维离散小波可W提取眼底图像的低频系数、水平高频系数、垂直高频系数，对角高 频系数，但在传统的离散小波变换在每次计算小波系数的过程中会丢弃前一步的尺度和小 波结果，形成类金字塔样分解。而UWT在计算DWT过程中，使用一个NXN环形平移的矩阵 TX来代替原来的X矩阵。如果X矩阵是一维的矢量X=技。;Xi;... ;XN_2;XN_i]，则
[0053] TX=Rn_i;X〇;Xi;…；XnJ，T-iX=化;乂2;…；Xn-2;XJ
[0054] 然后DWT计算一次X，再计算一次TX，并将各自得到的DWT系数进行融合来形 成UWT的小波系数【D.B.Percival,A.T.Walden.WaveletMethodsforTimeSeries Analysis.CambridgeUniversityPress,Cambridge(2000)】。然后，通过对小波系数的重 构得到不同尺度下的视网膜图像。
[00巧]本实验中，通过UWT得到了视网膜图像不同尺度的小波系数，然后根据不同的目 的对小波系数进行相应的处理，然后再对处理后的小波系数进行图像重构从而完成图像的 处理。重构后的眼底图像就得到了增强。视网膜血管与视网膜图像背景的对比度明显增强， 图像清晰（图1)。
[0056] 本实验选择第四层的UWT重构后图像进行下一步的特征提取。
[0057] 本方法采用基于模糊C-均值聚类分类器的视网膜血管像素分类。模糊聚类分析 就是把一个模糊的不确定的没有类别标记的样本按照某种准则划分为若干个子群集，让相 似的样本尽可能归为一类，而把不相似的样本划分到不同的类中。是一种无监督模式识别。 模糊聚类算法是一种基于函数最优方法的聚类算法，需要使用微积分计算技术求最优代价 函数，在基于概率算法的聚类方法中使用概率密度函数。
[0058] 模糊集合最根本的方法是其所包含元素及其相应的隶属度用函数表示出来，模糊 集合的数学定义如下：
[0059] 称A是给定论域U上的一个模糊子集的含义是；对任何UGU都对应一个唯一的 实数yA(u)G[0，1]，有序对集合{u，yA(u)|uGU}均可作为U上的模糊集，yA(u)表示 U对A的隶属度。该表明构造了一个映射；
[0060]yA(u) :U^ [0, 1],u^yA(u)
[0061] 该个映射称为模糊集A的隶属度函数。一个定义在空间X= {x}上的隶属度函数 就定义了一个模糊集合A，或者叫定义在论域X= {x}上的模糊子集对于有限个对象XI， X2,……，x。模糊集合^可W表示为：
[0062] ^ = !(片,臣乂! ~〇
[0063] 模糊C-均值聚类（FCM)的核屯、思想是寻求合适的隶属度和聚类中屯、，使得聚类内 代价函数的方差和迭代误差最小，代价函数的值是像素到聚类中屯、L2测度的加权累积和， 表达形式为：
[0064]
[0065] J扣，V)表示区域的像素到聚类中屯、加权距离的平方和。Uik表示X i属于第k类区 域的隶属度。m为模糊指数，me[1，-]，表示控制分类矩阵U的模糊程度，是隶属度的加 权指数，决定分类结果的模糊程度，其值越大分类越模糊。同时满足约束条件：=1， z=l 求解此类极小值的问题，可W通过拉格朗日乘数法来求解。首先，由于矩阵U中各列独立， 则有：
[0066]
[0067] 求解极小值问题则转化为求解
的极小值，且在约束条件
.通过求解偏导数
可得出，在约束条件下，其隶属函数和聚类中屯、的更新依赖于W下两个函 数：
[ 0068]
[0069] 由该两个公式得出最佳模糊分类矩阵和聚类中屯、，该一过程可W通过迭代方法实 现[抽。
[0070] 本实验过程中，设定FCM的类别数为2,隶属度的加权指数m取2。最大迭代次数 为100次，判别误差为le-5,提取眼底图像中属于血管的像素点。
[0071] 本发明的方法具有较高的灵敏度化82)，算法运行迭代至函数收敛次数平均为 50次左右，实现分割的平均时间13. 4s。
[0072] 步骤二，视网膜血管的骨架化和拓扑分层
[007引如图2所示，本发明首先按照图像4邻域法则，识别出眼底血管网络上的特征点如 分叉点和终末点。
[0074] 然后根据建立的特征点，得到两个特征点之间的曲线部分，并根据图2的层次标 定方法，对眼底血管骨架网络进行重新标定，标定的结果如图3所示。
[00巧]本方法提出了一种定义血管树的分支的编号方法，该理论定义所有的子结点都是 由根节点派生出来的。结点是按由上往下，由根结点到子节点编号的，首先将根编号发出的 分支编号为1 ;然后根节点向下发出的子节点的分支编号+1，W此类推。（图2)
[0076] 步骤S，视网膜分层的迂曲度计算
[0077] 本发明根据血管曲线的每个像素点的X，Y坐标作为输入值，利用协方差矩阵的迹 和曲线的像素点个数对血管迂曲度进行建模。
[0078] 在本发明中，定义眼底血管待测曲线的每个像素点X，Y坐标为两个随机变量X和 Y，根据协方差（记为Cov狂，Y))定义：
[0079]Cov化Y) =E{技-E狂)][Y-E(Y) ]}
[0080]其中，E佩二化+X2+. . .Xn)/n;E(Y)=化+Y2+. . . +YJ/n
[0081]
陈为rv〇(，Y)的协方差矩阵，该里X，Y分别对应于待测血管曲线上的每 个像素的横纵坐标值。
[0082] 本发明中定义迂曲度化r:
[0083]Tur=n]/(C11+C22)
[0084] 该式中n为待测曲线上的像素点个数，即骨架化后图中像素点个数n，C11+C22为该 图骨架化后所有像素点X，Y坐标值的协方差矩阵的迹。
[0085] 为了验证我们的算法，我们构建出不同的迂曲模型，如图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 峰两谷系列中，第S峰较高，测得迂曲度为10. 8641，S峰齐 高，测得迂曲度为11. 2867,第一峰较高，测得迂曲度为11. 1562。
[0086]如图5所示，取一近半圆曲线的原图，将其分别旋转0°、15°、90°，分别计算 其迂曲度，S者均为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。
[0087]W上结果说明本发明的迂曲度算法能够较好地表征出模型中的迂曲情况。
[0088] 实施例2;运用本发明方法测量早期糖尿病患者视网膜血管改变
[0089]W10幅正常视网膜图像和10幅糖尿病视网膜图像为实验对象，图像有临床眼底 镜工作站（南通大学附属医院）提供。图像来自于20例临床确诊为2型糖尿病患者，诊断 标准为《1997年美国糖尿病协会提出的糖尿病诊断标准》1、有糖尿病症状，并且任意血糖 ^ 11.Immol/L。2、空腹血糖 ^ 7.Ommol/L。3、糖耐量试验 2 小时血糖 ^ 11.Immol/L。0GTT 仍按世界卫生组织的要求进行。由于有10例合并高血压，其眼底照片被排除。最终纳入10 幅图像，其中男性6例，女性4例，平均年龄为51. 2岁，10例眼底图像按照国际临床DR的 严重程度分级和我国糖尿病视网膜病变的分类方法，其中7例被诊断为无明显糖尿病视网 膜病变值i油eticretinopathy,DR)，另外3例被诊断出轻度非增生性DR(有微血管瘤，单 纯型I期）。从每一幅图像中选择一段血管枝进行新迂曲度计算。先得到视网膜血管骨架 化后的矩阵，然后利用本发明的算法可W实现自动分层标记。在实际血管树图像处理中，会 经常发现出血管现交叉情况，对于该类图像，我们设计了一种半自动的血管分层编程方法。 首先将眼底血管图像进行分解，其中每个部分内部只有分叉类型，另外图中保留起始端的 部分血管枝，便于最后重新组合。再对分解的图进行组合，将其中血管支数据整合成一组数 据，再进行使用。
[0090] 如图6所示，从折线图中可w看出糖尿病眼底血管支的迂曲度起伏不定，最后趋 于增大，而正常血管支的迂曲度更加平稳，波动较小。
[0091] 可见，本发明方法发现早期糖尿病患者在高的视网膜层次上迂曲度提高。
[0092] W上已对本发明创造的较佳实施例进行了具体说明，但本发明创造并不限于所述 实施例，熟悉本领域的技术人员在不违背本发明创造精神的前提下还可作出种种的等同的 变型或替换，该些等同的变型或替换均包含在本申请权利要求所限定的范围内。
【主权项】
1. 一种基于眼底镜图像的视网膜血管迂曲度计算方法，包括以下步骤： 步骤一：首先使用数字化眼底镜获取筛查人群的眼底图像； 步骤二：利用非下采样的离散小波变换对图像进行增强； 步骤三：然后视网膜灰度图像局部熵的纹理提取； 步骤四：然后利用模糊C聚类的方法对视网膜血管进行分割； 步骤五：对分割后的血管进行骨架化； 步骤六：计算骨架的拓扑层次； 步骤七：利用以下迂曲度计算模型对血管骨架每个层次进行迂曲度计算： 根据血管曲线的每个像素点的X，Y坐标作为输入值，利用协方差矩阵的迹和曲线的像 素点个数对血管迂曲度进行建模； 眼底血管待测曲线的每个像素点X，Y坐标为两个随机变量X和Y，根据协方差记为 Cov (X，Y)定义： Cov (X, Y) = E {[X-E (X) ] [Y-E (Y) ]} 其中，E(X) = (XJX2+. · · Xn)/n ;E(Y) = (YJY2+. · · +Yn)/n称为rv(X，Y)的协方差矩阵，这里X，Y分别对应于待测血管曲线上的每个像 素的横纵坐标值； 迂曲度Tur定义： Tur = n2/ (Cn+C22) 该式中η为待测曲线上的像素点个数，即骨架化后图中像素点个数n，Cn+C2#该图骨 架化后所有像素点X，Y坐标值的协方差矩阵的迹。2. -种如权利要求1所述的基于眼底镜图像的视网膜血管迂曲度计算方法在糖尿病 早期计算机筛查模型或系统中的应用。
【专利摘要】本发明属于医学图像处理及应用领域，提供了一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用。本发明首先使用数字化眼底镜获取筛查人群的眼底图像，再利用非下采样的离散小波变换(UDWT)对图像进行增强；然后视网膜灰度图像局部熵的纹理提取，利用模糊C聚类(FCM)的方法对视网膜血管进行分割；最后对分割后的血管进行骨架化，并计算骨架的拓扑层次，并用本发明的迂曲度计算模型对血管骨架进行迂曲度计算。本发明的方法实施简单，精度可靠，便于临床应用。
【IPC分类】G06T7/00
【公开号】CN104881862
【申请号】CN201510159636
【发明人】吴辉群, 董建成, 耿兴云, 张晓峰
【申请人】南通大学
【公开日】2015年9月2日
【申请日】2015年4月3日