眼部疾病是危害人类健康的一类重大疾病。常见的眼部疾病包括青光眼、糖尿病性视网膜病变、黄斑病变等(Thakur和Juneja，2018)。其中青光眼疾病是致盲的主要原因之一, 而且青光眼引起的视觉损伤是不可逆的, 后果极为严重。早期及时发现治疗可极大地降低视力损伤, 使多数患者终生保持有用的视功能(Almazroa等，2015)。青光眼的早期检查通常是通过彩色眼底图像中视杯(optic cup, DC)垂直直径和视盘(optic disc, OD)垂直直径的长度比例即杯盘比(cup to disc ratio, CDR)关系确定。图 1展示了眼底图像视杯、视盘及其比例关系。CDR已经成为当今青光眼疾病诊断和患病程度的一个重要衡量指标。因此, 在眼底图像中, 如何准确地分割视杯和视盘区域进而准确地计算CDR获得了广泛关注(Hagiwara等，2018；Jonas等，2000)。传统的视杯和视盘区域分割通常是由有经验的医生手动完成, 这种分割方式耗费大量精力并且分割效率较低。随着图像处理和计算机视觉技术的发展, 如何通过图像处理技术，实现眼底图像中自动准确地分割视杯和视盘区域已经成为了一个重要课题。

视杯和视盘的分割方法有基于超像素(Cheng等，2013)、血管空间结构先验分析(Joshi等，2011；Guo等，2019)、聚类分析(Khalid等，2014)和活动轮廓模型(Zhou等，2010)等。Cheng 等人(2013)提出了基于超像素的视杯视盘分割方法，在视盘的分割中, 采用直方图和中心环绕统计方法判断每个超像素分割结果为视盘区域或非视盘区域，在此基础上，通过先验位置信息实现视杯区域分割。Joshi 等人(2011)考虑视杯附近血管的弯曲特性，通过获得血管弯曲的可靠集，拟合期望的杯盘边界。Guo 等人(2019)首先应用一组面向血管几何结构的COSFIRE(combination of shifted filter responses)滤波器定位血管的主要分叉点，在此基础上获得包含视杯和视盘的感兴趣区域，应用一组面向视盘几何结构的COSFIRE滤波器检测视盘区域。最后，应用基于广义矩阵学习向量量化的分类方法将视盘区域分割成视杯区域和盘沿区域，将视杯和视盘的边界近似成椭圆形状。Khalid 等人(2014)利用改进模糊C均值方法分割眼底图像中的视杯和视盘区域。Noor 等人(2013)首先提取视盘附近感兴趣区域，在感兴趣区域内通过对彩色眼底图像的多通道设置不同阈值进行视杯视盘的分割。Zheng 等人(2013)构建了一个基于图割法的全局优化框架，结合视杯和视盘的形状、位置和几何关系等多种先验信息，提出视杯视盘分割问题的可靠解决方案。Chakravarty和Sivaswamy(2017)提出一种基于边界条件随机场方法，可有效模拟视杯和视盘边界之间的空间结构关系，同时实现视杯和视盘边界的提取。

活动轮廓模型广泛应用于眼底图像中视杯和视盘的轮廓提取，获得了很好的效果。Zhou 等人(2010)改进传统的梯度向量流活动轮廓模型(gradient vector field-Snake, GVF-Snake)以实现视盘边界的准确分割，在标准GVF力场中引入均值漂移算子, 使得活动轮廓在内力和外力的共同作用下向视盘边界演化，准确提取视盘边界。Esmaeili 等人(2012)采用傅里叶相关系数滤波法获得视盘初始边界, 以此作为正则化水平集演化方法(distance regularized level set evolution, DRLSE)的初始轮廓, 分割视盘边界, 获得准确结果。Dai 等人(2017)提出了一种基于形状约束与区域约束拟合的变分活动轮廓模型用于自动分割眼底图像中的视盘区域。Gao 等人(2018)在能量模型中引入视盘的先验椭圆形状约束，提高分割算法的精度。在眼底图像的视杯分割中, 由于视杯边界附近存在大量的血管结构干扰, 使得视杯区域的分割比视盘区域更具挑战性。此外, 神经视网膜区域的灰度是不均匀的，进一步增加了视杯边界提取的复杂性。Wong 等人(2008)提出一种基于边缘的变分水平集方法提取眼底图像视盘边界，根据视盘边界提取结果，采用椭圆拟合进一步提取视杯边界。Mittapalli和Kande(2016)采用局部二值拟合模型提取视盘轮廓，在此基础上，依据视杯的结构和灰度特性分割视杯区域。Lotankar 等人(2015)通过测地线活动轮廓模型对视盘区域进行分割，接着使用CMY(cyan, magenta, yellow)色彩空间中M(magenta)通道的颜色特征信息对视杯区域进行提取。郑姗等人(2014)提出了一种基于先验形状约束的多相Chan-Vese(C-V)模型, 该方法使用两个水平集函数分割视盘区域和视杯区域, 并且约束水平集函数在演化过程中始终保持一定的形状。该方法能够实现视杯和视盘区域同时分割, 并取得了较理想的分割结果。但是该方法对灰度不均匀图像处理效果不理想, 并且同时使用两个水平集函数导致轮廓演化速度较慢。Thakur和Juneja(2019)提出一种基于聚类融合的水平集方法应用于视杯和视盘的分割。该方法融合基于自适应正则化核的模糊C均值聚类和直觉模糊C均值聚类两种聚类方法提取视杯和视盘区域，在此基础上以分类结果作为初始轮廓，采用水平集方法分别对视杯和视盘边界进行进一步修正，获得准确的视杯和视盘边界。

上述大部分算法虽然成功实现了眼底图像视盘及视杯区域的分割，但通常没有考虑视杯和视盘的解剖几何结构关系，且算法易受视网膜病理区域、血管覆盖和图像的灰度不均匀性等因素的影响。针对眼底图像的视盘及视杯分割问题，本文提出一种基于双层水平集描述的眼底图像视杯视盘分割方法。该方法可以通过一个水平集函数的2个不同层级分别代表视网膜眼底图像中的视杯和视盘边界。

本文方法主要有3方面贡献：1)考虑了视杯与视盘间的距离关系，应用一个水平集函数可实现具有不同区域亮度的视杯和视盘区域的分割；2)针对图像灰度不均匀问题进行改进；3)降低血管结构对视盘视杯分割结果的影响。

本文通过视杯与视盘分割时水平集函数的整体演化过程分析，验证了采用一个水平集函数实现视杯和视盘区域准确分割的有效性。通过与基于COSFIRE滤波模型的视杯视盘分割方法、基于先验形状约束的多相C-V模型和基于聚类融合的水平集方法的对比实验，验证本文方法能够有效克服眼底图像灰度不均匀性的影响和血管及病变区域的干扰，更加准确地分割视杯和视盘区域。

在视网膜眼底图像中，准确获取视杯和视盘的视网膜结构，是实现青光眼诊断的前提。通过对视网膜眼底图像及其视杯和视盘区域进行大量实验和观察，可以发现视网膜眼底图像有如下特点：1)在灰度图像中，视网膜眼底图像是灰度不均匀的，并且视盘区域和视杯区域具有不同的区域亮度；2)视杯和视盘的轮廓近似两个椭圆，视杯位于视盘内部，视盘与视杯轮廓之间的距离变化是均匀的；3)视盘的轮廓特征较为明显，由于视盘中包含许多病灶区域，并且视杯边界附近存在大量的血管结构，使得视杯的轮廓存在不连续现象。

由于视盘和视杯区域具有不同的区域亮度, 因此现有大部分方法都采用2个水平集函数对视杯和视盘轮廓进行提取，导致算法的计算量较大，迭代收敛速度较慢。视杯区域完全位于视盘内部，视盘与视杯轮廓不相交，并且彼此之间的距离变化均匀。因此，可以采用一个水平集函数的2个不同层分别描述视杯和视盘轮廓。本文借鉴Feng等人(2016)提出的分层描述不同轮廓的思想，采用一个水平集函数的0层和$k$层分别描述视杯和视盘轮廓。假设视网膜眼底图像${\mathit{\boldsymbol{I}}}$在定义域${\mathit{\boldsymbol{Ω}}}$上是连续变化的, $ϕ$表示水平集函数，$C_{0}$和$C_{k}$分别代表水平集函数的第0层和第$k$层。因此，该图像在其定义域${\mathit{\boldsymbol{Ω}}}$上被分为3部分，定义为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_1} \buildrel \Delta \over = \{ x:\phi (x) < 0\} }\\ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_2} \buildrel \Delta \over = \{ x:0 < \phi (x) < k\} }\\ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_3} \buildrel \Delta \over = \{ x:\phi (x) > k\} } \end{array}} \right. $

(1)

引入Heaviside函数(Zhao等，1996)，获得各子区域，即

$ \left\{ {\begin{array}{*{20}{l}} {{M_1}(\phi (x)) = 1 - H(\phi (x))}\\ {{M_2}(\phi (x)) = H(\phi (x)) - H(\phi (x) - k)}\\ {{M_3}(\phi (x)) = H(\phi (x) - k)} \end{array}} \right. $

(2)

式中, ${\mathit{\boldsymbol{M}}}_{1}$, ${\mathit{\boldsymbol{M}}}_{2}$, ${\mathit{\boldsymbol{M}}}_{3}$分别表示视杯，视盘和背景区域。$H$表示Heaviside函数。

为保证水平集函数的0层和$k$层之间距离缓慢平滑变化，需要构建一个距离约束项进行约束，为

$ \begin{array}{l} D(\phi, \alpha) = {\mu _1}\int {\frac{1}{2}} {(|\nabla \phi (x)| - \alpha (x))^2}{\rm{d}}x + \\ \ \ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mu _2}\int {\frac{1}{2}} {(\nabla \alpha (x))^2}{\rm{d}}x \end{array} $

(3)

式中，Δ表示梯度算子，$μ_{1}>0$，$μ_{2}>0$为权重系数，$α(x)$为一个光滑函数，通常情况下$α(x)$的初始值为1。第1项约束水平集函数的梯度为一个光滑函数，逼近函数$α(x)$，第2项保证函数$α(x)$的光滑性。

为保证水平集函数的0层和$k$层在演化过程中时刻保持光滑，本文定义长度项为

$ \begin{array}{*{20}{c}} {L(\phi) = {v_1}\int | \nabla H(\phi (x))|{\rm{d}}x + }\\ {{v_2}\int | \nabla H(\phi (x) - k)|{\rm{d}}x} \end{array} $

(4)

式中，第1项和第2项分别为水平集函数0层和$k$层的长度约束项，$v_{1}$和$v_{2}$分别为第1项和第2项的权重系数。

视网膜眼底图像通常是灰度不均匀的，需要构建数据驱动项克服图像灰度不均匀带来的影响, 进而准确提取视杯和视盘轮廓。针对图像灰度不均匀问题，一些充分考虑图像局部信息的活动轮廓模型如局部二值拟合模型(local binary fitting，LBF)(Li等，2007)、局部强度聚类模型(local intensity clustering，LIC)(Li等，2011)等方法相继提出, 取得了满意的结果。本文采用LIC模型的数据驱动项形式，定义为

$ \begin{array}{*{20}{c}} {E(\phi, {c_1}, {c_2}, {c_3}, b(x)) = }\\ {\left. {\int {\left({\sum\limits_{i = 1}^3 {{\lambda _i}} } \right.\int {{K_\rho }} (x - y)|I(x) - b(y){c_i}{|^2}{M_i}(} \phi (x)){\rm{d}}x} \right){\rm{d}}y} \end{array} $

(5)

式中，$b$表示灰度偏移场，$c_{i}$为一个分片常量函数，$bc_{i}$拟合图像的局部灰度加权均值, $λ_{i}$为权重参数。函数$K_{ρ}$是截断核函数, 定义为

$ {K_\rho }(u) = \left\{ {\begin{array}{*{20}{l}} a&{|u| \le \rho }\\ 0&{|u| > \rho } \end{array}} \right. $

(6)

式中，$u$表示距离值，$ρ$表示核函数半径，$a>0$以满足$\int_{|u| \le \rho } {{K_\rho }} (u) = 1$。

在眼底图像中，视杯位于视盘内部，轮廓均近似两个椭圆形状。视盘中包含许多病灶区域，并且视杯边界附近存在大量的血管结构，这对轮廓精确提取造成极大干扰。因此，需要根据视杯和视盘的形状特点，构建一个先验形状约束项，通过引入视杯和视盘的椭圆形状先验约束，使得水平集函数在演化过程中，0层和$k$层轮廓总是逼近椭圆形状，有效克服血管和病灶区域对轮廓提取的影响。

视杯轮廓椭圆形状先验约束项定义为

$ \begin{array}{*{20}{c}} {{\phi _1}(x) = 1 - \left[ {\frac{{{{((x - {x_1}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1} + (y - {y_1}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1})}^2}}}{{a_1^2}}} \right. + }\\ {{{\left. {\frac{{{{(- (x - {x_1}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1} + (y - {y_1}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1})}^2}}}{{b_1^2}}} \right]}^{\frac{1}{2}}}} \end{array} $

(7)

视盘轮廓椭圆形状先验约束项定义为

$ \begin{array}{*{20}{c}} {{\phi _2}(x) = 1 - \left[ {\frac{{{{((x - {x_2}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2} + (y - {y_2}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2})}^2}}}{{a_1^2}}} \right. + }\\ {{{\left. {\frac{{{{(- (x - {x_2}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2} + (y - {y_2}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2})}^2}}}{{b_2^2}}} \right]}^{\frac{1}{2}}}} \end{array} $

(8)

式中，$(x_{1}, y_{1})$和$(x_{2}, y_{2})$分别为两个椭圆中心的坐标，$(a_{1}, b_{1})$和$(a_{2}, b_{2})$分别为两个椭圆的半长轴(或半短轴)和半短轴(或半长轴)，$θ_{1}$和$θ_{2}$为旋转角度。本文中将椭圆形状先验约束中涉及到的参数表示为$(a, b, θ, x, y)$。$ϕ_{1}(x)$，$ϕ_{2}(x)$是用于描述椭圆形状的水平集函数。本文定义如下的先验形状约束项，即

$ \begin{array}{l} {E_{{\rm{ prior }}}} = {\gamma _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {(H(} \phi (x)) - H{({\phi _1}(x))^2}{\rm{d}}x + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\gamma _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {(H(} \phi (x) - k) - H{({\phi _2}(x))^2}{\rm{d}}x \end{array} $

(9)

式中，第1项为视杯形状约束项，第2项为视盘形状约束项。$γ_{1}>0$，$γ_{2}>0$为权重系数。

基于上述分析，本文算法能量泛函形式为

$ {E_{{\rm{ total }}}}(\phi) = D(\phi) + L(\phi) + E(\phi) + {E_{{\rm{ prior }}}}(\phi) $

(10)

式中，$D$是距离正则项，能够保证水平集函数的0层和$k$层之间的距离变化是光滑的；$L$是长度正则项，使得水平集函数的0层和$k$层在演化过程中时刻保持光滑。$E$是数据驱动项，源于图像信息，使得模型能够克服图像灰度不均匀的影响。$E_{\rm prior}$是先验形状约束项，使得水平集函数在演化过程中, 0层和$k$层轮廓总是逼近椭圆形状，可有效克服血管和病灶区域对轮廓提取的影响。

眼底图像中视杯及视盘轮廓提取是通过最小化总能量$E_{\rm total}$实现的。在$c_{1}$, $c_{2}$, $c_{3}$, $α$, $b(x)$, $a$, $b$, $θ$, $x$和$y$固定的条件下, 采用梯度下降法, 水平集函数$ϕ$的演化方程为

$ \begin{array}{l} \begin{array}{*{20}{c}} {\frac{{\partial \phi }}{{\partial t}} = {\lambda _1}\delta (\phi){e_1} - {\lambda _2}(\delta (\phi) - \delta (\phi - k)){e_2} - }\\ {{\lambda _3}\delta (\phi - k){e_3} + } \end{array}\\ \begin{array}{*{20}{c}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({v_1}\delta (\phi) + {v_2}\delta (\phi - k)) {\rm{div}} \left({\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mu \left({{\nabla ^2}\phi - \alpha {\rm{div}} \left({\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right)} \right) - } \end{array}\\ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{\gamma _1}\delta (\phi)(H(\phi (x)) - H({\phi _1}(x))) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{\gamma _2}\delta (\phi)(H(\phi (x)) - H({\phi _2}(x)))} \end{array} \end{array} $

(11)

式中，$δ(x)$函数是Heaviside函数的导数，div(·)为散度算子，$e$为

$ {e_i}(x) = \int {{K_\rho }} (x - y)|I(x) - b(y){c_i}{|^2}{\rm{d}}y $

的缩写。

在$ϕ$, $α$, $b(x)$, $a$, $b$, $θ$, $x$和$y$固定的条件下, 最小化能量泛函$E_{\rm total}$，$c_{1}$，$c_{2}$，$c_{3}$可分别定义为

$ \left\{ \begin{array}{l} {c_1} = \frac{{\int {(b(} x) * {K_\rho })\mathit{\boldsymbol{I}}{M_1}(\phi (y)){\rm{d}}y}}{{\int {(b(} x{)^2} * {K_\rho }){M_1}(\phi (y)){\rm{d}}y}}\\ {c_2} = \frac{{\int {(b(} x) * {K_\rho })\mathit{\boldsymbol{I}}{M_2}(\phi (y)){\rm{d}}y}}{{\int {(b(} x{)^2} * {K_\rho }){M_2}(\phi (y)){\rm{d}}y}}\\ {c_3} = \frac{{\int {(b(} x) * {K_\rho })\mathit{\boldsymbol{I}}{M_3}(\phi (y)){\rm{d}}y}}{{\int {(b(} x{)^2} * {K_\rho }){M_3}(\phi (y)){\rm{d}}y}} \end{array} \right. $

(12)

在$ϕ$，$c_{1}$，$c_{2}$，$c_{3}$，$b(x)$，$a$，$b$，$θ$，$x$和$y$固定的条件下，采用梯度下降法，最小化距离正则项$D(ϕ, α)$，得到$α$的演化方程为

$ \frac{{\partial \alpha (x)}}{{\partial t}} = \mu (|\nabla \phi (x)| - \alpha (x)) + 2\omega {\kern 1pt} {\kern 1pt} {\kern 1pt} {\nabla ^2}\alpha (x) $

(13)

在$ϕ$，$α$，$c_{1}$，$c_{2}$，$c_{3}$，$a$，$b$，$θ$，$x$和$y$固定的条件下，最小化能量泛函$E_{\rm total}$，$b(x)$的演化方程为

$ b(x) = \frac{{\sum\limits_{i = 1}^3 {{K_\rho }} * (\mathit{\boldsymbol{I}}{M_i}(\phi (x))){c_i}}}{{\sum\limits_{i = 1}^3 {{K_\rho }} * {M_i}(\phi (x))c_i^2}} $

(14)

在$ϕ$，$c_{1}$，$c_{2}$，$c_{3}$，$α$和$b(x)$固定的条件下，最小化能量泛函$E_{\rm total}$，$a, b, θ, x, y$的演化方程为

$ {\frac{{{\rm{d}}{a_1}}}{{{\rm{d}}t}} = 2{\gamma _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 1}}} \left({\frac{{A_1^2}}{{a_1^3}}} \right){\rm{d}}x} $

(15)

$ {\frac{{{\rm{d}}{b_1}}}{{{\rm{d}}t}} = 2{\gamma _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 1}}} \left({\frac{{B_1^2}}{{b_1^3}}} \right){\rm{d}}x} $

(16)

$ \frac{{{\rm{d}}{x_1}}}{{{\rm{d}}t}} = 2{\gamma _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 1}}} \left({\frac{{{A_1}{\kern 1pt} {\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1}}}{{a_1^2}} - \frac{{{B_1}{\kern 1pt} {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1}}}{{b_1^2}}} \right){\rm{d}}x $

(17)

$ \frac{{{\rm{d}}{y_1}}}{{{\rm{d}}t}} = 2{\gamma _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 1}}} \left({\frac{{{A_1}{\kern 1pt} {\kern 1pt} {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1}}}{{a_1^2}} - \frac{{{B_1}{\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1}}}{{b_1^2}}} \right){\rm{d}}x $

(18)

$ \frac{{{\rm{d}}{\theta _1}}}{{{\rm{d}}t}} = 2{\gamma _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 1}}} \left({\frac{{ - {A_1}{\kern 1pt} {\kern 1pt} {B_1}}}{{a_1^2}} + \frac{{{A_1}{\kern 1pt} {\kern 1pt} {B_1}}}{{b_1^2}}} \right){\rm{d}}x $

(19)

$ \frac{{{\rm{d}}{a_2}}}{{{\rm{d}}t}} = 2{\gamma _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 2}}} \left({\frac{{A_2^2}}{{a_2^3}}} \right){\rm{d}}x $

(20)

$ \frac{{{\rm{d}}{b_2}}}{{{\rm{d}}t}} = 2{\gamma _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 2}}} \left({\frac{{B_2^2}}{{b_2^3}}} \right){\rm{d}}x $

(21)

$ \frac{{{\rm{d}}{x_2}}}{{{\rm{d}}t}} = 2{\gamma _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 2}}} \left({\frac{{{A_2}{\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2}}}{{a_2^2}} - \frac{{{B_2}{\kern 1pt} {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2}}}{{b_2^2}}} \right){\rm{d}}x $

(22)

$ \frac{{{\rm{d}}{y_2}}}{{{\rm{d}}t}} = 2{\gamma _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 2}}} \left({\frac{{{A_2}{\kern 1pt} {\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2}}}{{a_2^2}} - \frac{{{B_2}{\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2}}}{{b_2^2}}} \right){\rm{d}}x $

(23)

$ \frac{{{\rm{d}}{\theta _2}}}{{{\rm{d}}t}} = 2{\gamma _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{H_{\ell 2}}} \left({\frac{{ - {A_2}{\kern 1pt} {B_2}}}{{a_2^2}} + \frac{{{A_2}{\kern 1pt} {B_2}}}{{b_2^2}}} \right){\rm{d}}x $

(24)

式中，$H_{\ell 1}$，$H_{\ell 2}$，$A_{1}$，$A_{2}$，$B_{1}$，$B_{2}$定义为

$ \begin{array}{*{20}{c}} {{H_{\ell 1}} = (H(\phi) - H({\phi _1}))\delta ({\phi _1})}\\ {{H_{\ell 2}} = (H(\phi) - H({\phi _2}))\delta ({\phi _2})}\\ {{A_1} = (x - {x_1}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1} + (y - {y_1}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1}}\\ {{B_1} = - (x - {x_1}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1} + (y - {y_1}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _1}}\\ {{A_2} = (x - {x_2}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2} + (y - {y_2}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2}}\\ {{B_2} = - (x - {x_2}){\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2} + (y - {y_2}){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _2}} \end{array} $

在实际计算中，为保持数值计算的稳定性，使用正则化Heaviside函数$H_{ε}$，定义为

$ {H_\varepsilon } = \frac{1}{2}\left[ {1 + \frac{2}{\pi }{\rm{arctan}}\left({\frac{x}{\varepsilon }} \right)} \right] $

(25)

$H_{ε}$的导数Dirac函数$δ_{ε}(x)$定义为

$ {\delta _\varepsilon }(x) = \frac{1}{\pi }\frac{\varepsilon }{{{\varepsilon ^2} + {x^2}}} $

(26)

本文采用显式方案进行数值计算求解，式(11)可离散为

$ \begin{array}{l} \begin{array}{*{20}{c}} {{\phi ^{n + 1}} = {\phi ^n} + \Delta t({\lambda _1}\delta ({\phi ^n}){e_1} - {\lambda _2}(\delta ({\phi ^n}) - }\\ {\delta ({\phi ^n} - k)){e_2} - {\lambda _3}\delta ({\phi ^n} - k){e_3}) + ({v_1}\delta ({\phi ^n}) + } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{v_2}\delta ({\phi ^n} - k)) {\rm{div}} \left({\frac{{\nabla {\phi ^n}}}{{|\nabla {\phi ^n}|}}} \right) + }\\ {\mu \left({{\nabla ^2}{\phi ^n} - \alpha {\rm{div}} \left({\frac{{\nabla {\phi ^n}}}{{|\nabla {\phi ^n}|}}} \right)} \right) - } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {2{\gamma _1}\delta ({\phi ^n})(H({\phi ^n}(x)) - H(\phi _1^n(x))}\\ {2{\gamma _2}\delta ({\phi ^n})(H({\phi ^n}(x)) - H(\phi _2^n(x)))} \end{array} \end{array} $

(27)

式中，$Δt$为迭代时间步长，$n$为迭代次数，通过对上式迭代求解即可获得最终分割结果。

在水平集函数的演化过程中，本文采用交互扩散方法(Zhang等，2013)解决水平集函数重新初始化问题，即

$ {\phi ^{n + 1}} = {\phi ^n} + \Delta {t_1} \cdot \Delta {\phi ^n} $

(28)

本文算法具体步骤如下：

1) 初始化水平集函数$ϕ(x)$，$ϕ_{1}(x)$和$ϕ_{2}(x)$为符号距离函数；初始化函数$α(x)$；初始化$b(x)$。

2) 根据式(12)计算$c_{1}$，$c_{2}$，$c_{3}$。

3) 根据式(15)—式(24)计算更新$a$，$b$，$θ$，$x$，$y$。

4) 根据式(11)演化水平集。

5) 根据式(28)平滑水平集函数。

6) 根据式(13)计算$α(x)$。

7) 根据式(14)计算$b(x)$。

8) 判断水平集函数收敛条件，本文采用如下准则判断水平集函数是否收敛，定义为

$ \begin{array}{*{20}{l}} {Q = |L({C_0}(n)) - L({C_0}(n - 1))| + }\\ \ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} |{\kern 1pt} {\kern 1pt} L({C_k}(n)) - L({C_k}(n - 1))|} \end{array} $

(29)

式中，$Q$为相邻两次迭代中0层和$k$层长度改变量的总和。$C_{0}$表示0层的闭合曲线，$C$$_{k}$表示$k$层的闭合曲线。如果$Q<τ$，$τ$为阈值，则水平集函数停止演化，否则，返回步骤2)继续执行。

采用本文方法分割眼底图像视盘及视杯区域，并与算法1：基于COSFIRE滤波模型的视杯视盘分割方法(Guo等，2019)；算法2：基于先验形状约束的多相C-V模型(Zheng等，2014)；算法3：基于聚类融合的水平集方法(Thakur和Juneja，2019)的分割结果进行对比，验证算法的有效性。本文使用的眼底图像来自CDRISHTI-GS1数据集(Sivaswamy等，2015), 该数据集由印度Aravind眼科医院收集并注释，主要用来验证应用计算机辅助诊断算法分割视杯及视盘的鲁棒性和有效性。该数据集共包含101幅彩色眼底图像。每幅图像的视野范围为30°，分辨率为2 896×1 944像素。数据集提供了视杯和视盘边界的真实值，该值由4位具有不同临床治疗经验的医生人工标记轮廓计算平均值得到。

为了评估本文所提眼底图像视杯与视盘分割算法的准确性，采用相似性度量函数Dice系数(Thakur和Juneja，2018)计算真实目标区域${\mathit{\boldsymbol{X}}}$与分割算法所获得的目标区域${\mathit{\boldsymbol{Y}}}$的重叠率。Dice系数定义为

$ {D_i} = \frac{{2 \times |\mathit{\boldsymbol{X}} \cap \mathit{\boldsymbol{Y}}|}}{{|\mathit{\boldsymbol{X}}| + |\mathit{\boldsymbol{Y}}|}} $

(30)

式中，|· |表示区域内的像素个数。Dice系数的取值范围是0到1，该值越大代表算法的准确性越高。

采用F1分数指标评估分割算法的分割效果，即

$ F1 = 2\frac{{P \cdot R}}{{P + R}} $

(31)

F1分数取值范围同样是在0到1之间，其值越大则分割效果越好。式中，$P = \frac{{|\mathit{\boldsymbol{X}} \cap \mathit{\boldsymbol{Y}}|}}{{|\mathit{\boldsymbol{X}}|}} $表示分割算法的准确性，$R = \frac{{|\mathit{\boldsymbol{X}} \cap \mathit{\boldsymbol{Y}}|}}{{|\mathit{\boldsymbol{Y}}|}}$表示分割算法的灵敏性。

本文实验参数设置为$Δt$=0.1，$Δt_{1}=0.01$，$ρ=3$，$μ=10$，$ω=0.5$，$λ_{1}=0.005$，$λ_{2}=0.006$，$λ_{3}=0.004$，$γ_{1}=100$，$γ_{2}=100$，$v_{1}=0.17×255×255$，$v_{2}=0.12×255×255$，$α=1$，$b$=1。

图 2展示了本文方法对患有青光眼疾病眼底图像中的视杯与视盘区域进行分割的实验过程, 第1行为视杯与视盘轮廓曲线的演化迭代过程, 第2行为演化过程对应的水平集函数, 图中绿色曲线表示视杯轮廓, 红色曲线表示视盘轮廓。通过水平集函数的演化过程可以观察到, 水平集函数一直保持光滑, 具有数值计算的稳定性, 保证了最终分割结果的准确性和稳定性。

图 3展示了在青光眼疾病的不同发展阶段中本文方法、算法1、算法2和算法3的分割实验结果。图 3包含了青光眼疾病早期、中期和晚期的眼底图像。通过对上述实验结果分析可知，算法1能够分割出视杯和视盘的大部分区域，但是对视杯和视盘轮廓的提取不够精确。算法2由于没有考虑图像局部信息，致使分割结果受图像的灰度不均匀性影响较大。算法3没有考虑视杯和视盘的形状信息，因此分割结果受血管和病变的干扰较为严重, 并且存在过分割和欠分割问题。与上述方法相比, 本文方法可有效克服血管的干扰, 拟合出视杯和视盘的椭圆形状轮廓，同时保持了视杯与视盘间距离的缓慢平滑变化特性, 获得了更准确的分割结果。这是因为本文方法的能量泛函中增加了视杯与视盘轮廓的先验形状信息约束项，降低了血管对轮廓平滑性的干扰，同时能量泛函中的距离正则项保证了视杯与视盘的距离关系。

图 4展示了应用本文方法在分割图 2中眼底图像时获得的灰度偏差场$b(x)$和偏差校正的眼底图像。偏差校正的图像通过$I (x)/b(x)$计算获取。由图 4(c)可知, 在偏差校正的眼底图像中视杯、视盘和背景区域的灰度变得均匀且具有明显的可区分性。

图 5分别展示了上述原始图像和偏差校正眼底图像的灰度直方图分布, 通过对两幅灰度直方图的比较得出, 在偏差校正图像的直方图中存在3个分离良好的波峰, 分别对应眼底图像的背景、视杯和视盘区域。图像的灰度得到显著区分，从而实现了较好的分割效果。反之，由于偏差引起的灰度强度的重叠分布，原始图像的灰度直方图不具有这种良好的灰度峰值可分离特性。

为进一步验证本文方法的有效性，应用Dice系数作为视杯与视盘分割方法性能的定量评价标准，将本文方法与算法1、算法2和算法3进行对比，实验结果如表 1所示。本文算法分别获得了67.52 %的视杯平均重叠率，81.04 %的视盘平均重叠率，均高于其他3种算法。此外，还采用F1分数的平均值分析了上述4种方法的视杯和视盘分割结果，如表 2所示。由表 2可知，与其他3种方法相比，本文方法的视杯和视盘F1分数均为最高，并且视杯分割结果的F1分数标准差$σ_{\rm F\_OC}$和视盘分割结果的F1分数标准差$σ_{\rm F\_OD}$均为最小。上述实验进一步验证了本文方法在分割性能上优于其他3种方法，具有更高的分割精度。

本文提出一种基于双层水平集描述的眼底图像视杯视盘分割方法。该方法采用一个水平集函数的两个不同层级分别表示视杯边界和视盘边界，通过距离约束实现视杯与视盘间距离平滑变化，使用图像局部信息处理灰度不均匀，引入椭圆形状先验约束克服病灶和血管干扰，使得视杯和视盘的最终分割结果近似两个椭圆。实验证明，本文算法能够同时快速准确地分割视杯和视盘区域，取得较好的分割效果。

但是，本文方法仍有进一步改进的空间。在少数眼底图像中，图像的灰度不均匀较为严重，使得视杯与视盘的分割结果不理想。在之后的研究工作中，将研究如何将更多的图像信息如方差和局部熵等引入到数据项中，使之能够更加有效地克服灰度不均匀的影响，提高算法的准确性和鲁棒性，实现对严重灰度不均匀图像的精准分割。