CN105160660B - 基于多特征高斯拟合的活动轮廓血管提取方法及系统 - Google Patents
基于多特征高斯拟合的活动轮廓血管提取方法及系统 Download PDFInfo
- Publication number
- CN105160660B CN105160660B CN201510504854.8A CN201510504854A CN105160660B CN 105160660 B CN105160660 B CN 105160660B CN 201510504854 A CN201510504854 A CN 201510504854A CN 105160660 B CN105160660 B CN 105160660B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msubsup
- msup
- mfrac
- vessel
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20112—Image segmentation details
- G06T2207/20116—Active contour; Active surface; Snakes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20172—Image enhancement details
- G06T2207/20192—Edge enhancement; Edge preservation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30041—Eye; Retina; Ophthalmic
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
- Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于多特征高斯拟合的活动轮廓血管提取方法，同时还公开了一种相关的系统：使用局部相位血管增强算法对原始视网膜图像进行血管增强处理，以突出血管所在区域和降低图像灰度不均匀性；将血管增强结果图和原始灰度图中对应的像素值作为两个相互独立的随机变量，用于构建一种新的基于局部高斯拟合的能量泛函活动轮廓血管分割算法，同时引入一个正则项以惩罚血管轮廓的不平滑性和不合理的轮廓曲线的长度，并适当保持血管轮廓的局部细节；将该分割算法引入变分水平集求解框架中，实现血管轮廓的全自动提取。本发明能够较为准确全面地提取出视网膜图像中的血管组织，为眼部疾病的治疗提供有效的辅助。
Description
技术领域
本发明涉及一种图像处理领域的图像提取方法，具体涉及一种基于多特征高斯拟合的活动轮廓血管提取方法及系统。
背景技术
根据世界卫生组织的统计，各种血管类疾病已成为严重危害人类健康的三大疾病之一，其中便包括眼部视网膜血管相关疾病，且其具有较高的致残率。避免这种后果的常用方法是对患者进行早期的预防诊断，而准确诊断的关键是精确的血管提取。这不仅可以提高管腔狭窄、动脉瘤和血管钙化等血管疾病诊断结果的可靠性，对血管的介入治疗、手术计划制定和手术精确导航等临床应用也具有重要价值，而且可以为图像配准、三维重建、计算机辅助诊断等图像处理过程提供有用的信息。
由于眼部视网膜图像的背景信息较为复杂，且灰度不均匀性严重。因此，准确完整地提取视网膜图像中的各种不同粗细的血管组织是较为困难的具有挑战性的任务。直接对原始视网膜图像进行血管分割通常仅能提取相对较粗的血管，因而是一个较为粗略的结果，往往不能满足临床应用所需的精确度。因此，需要对原始图像进行血管增强处理，并将增强结果与原图像一块用于血管提取，通过多种图像特征的使用，能够取得自动分割更好的结果。
血管提取算法包括区域生长法、阈值法以及活动轮廓(active contour)等几类算法，其中活动轮廓算法能快速准确地获取子像素级别的分割结果而具有广阔的应用前景。活动轮廓算法根据分割过程中使用图像信息的不同，可分为基于边界(edge-based)的活动轮廓和基于区域(region-based)的活动轮廓。基于边界的算法使用图像梯度界定所需目标物体的边界，梯度的计算依赖于图像局部信息，进而导致基于边界的算法对图像噪声较为敏感。此外，求解过程中初始曲线位置偏离目标物体边界的严重程度也直接决定了分割结果的好坏，即距离目标物体边界越远，分割性能越差。这些不足导致该类算法只能分割具有较强灰度对比度的图像，而不能克服灰度不均匀性造成的问题。当前血管提取中，较为常用的是基于区域的活动轮廓模型，但它们仅仅考虑原始图像的灰度统计信息，并且较少涉及图像局部区域内的全面的灰度统计量(如灰度标准差)，使得分割结果不理想。
发明内容
为了解决上述技术问题，本发明提供了一种基于多特征高斯拟合的活动轮廓血管提取方法及系统，以解决视网膜图像中复杂背景信息和严重灰度不均匀特性导致的分割结果粗略，准确性低以及分割不完整等问题，从而准确完整地提取视网膜图像中具有不同粗细的各种血管组织。
为了达到上述目的，本发明提供了以下技术方案，其归纳总结为：在现有基于区域的活动轮廓模型中，引入血管增强算法，并将增强后的血管图像作为一种新的图像特征，用于血管提取算法中；具体而言如下：
本发明的基于多特征高斯拟合的活动轮廓血管提取方法包括：
通过血管增强滤波模块对视网膜图像进行血管增强处理，即采用局部相位(localphase)血管增强算法将视网膜图像中血管对应的区域突显出来，降低复杂背景信息和像素灰度不均匀性对分割结果造成的干扰；其中，增强后的血管图像(vesselness map)可以作为一种新的特征图像用于血管提取，即在血管图像中，像素灰度值越高，其代表血管的可能性越大；反之则越小为图像背景的可能性越高大；
通过局部高斯拟合能量泛函模块将血管图像和原始图像中的像素值作为两个相互独立的随机变量，构建一种基于二维高斯概率分布拟合的能量泛函活动轮廓模型，即在血管图像和原始视网膜图像中的像素值通过二维高斯分布统计模型模拟一定图像范围内的血管和背景区域内的图像统计特性，从而构建一个区域的活动轮廓模型。仅使用该能量泛函可能存在一定的不足，因此还需要通过一个正则约束项，惩罚血管轮廓曲线的不平滑性、不当的曲线长度、以及保持轮廓的局部细节；
通过变分求解框架模块使用变分水平集方法对能量泛函进行数学计算，即血管轮廓以隐式水平集的形式表示，然后通过梯度下降流和欧拉-拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的血管轮廓。
进一步地，上述的方法中，局部相位血管增强算法包括以下实施中的相关参数：设置中心频率为5π/7，带宽为2倍频(octave)，滤波范围为15×15，滤波方向为0°、45°、90°和135°，图像尺度系数为3，整合权重为3，正则系数为3。通过局部相位的血管增强可以显著突出血管所在的区域，同时减缓不相关背景信息的干扰。
进一步地，上述的方法中，局部高斯拟合能量泛函模块的能量泛函构建步骤如下：
在血管图像和原始视网膜图像中，对任一像素点局部邻域内的灰度概率分布通过二维高斯概率分布模型进行模拟：
P＝pi,x(I(y),V(y))＝pi,x(I(y))pi,x(V(y))；
其中，pi,x(I(y))为原始灰度图像中，中心像素点为x的邻域内，任一像素点y处的图像灰度I的概率分布，和为邻域内的灰度平均值和标准差；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
pi,x(V(y))，以及分别依次为增强后的血管图像V对应的概率分布、局部平均值和标准差；
将血管所在区域的像素和背景区域的像素区分开，因此血管区域的像素点中，像素的灰度值应尽可能地接近于灰度平均值，即灰度标准差尽可能趋近于0；转化为概率模型后为：
Πypi,x(I(y))pi,x(V(y))→max；
将上式的最大化求解转化为最小化，并引入基于距离的高斯权重 其中，σ为高斯函数中的标准差参数，同时考虑到增强后的血管图像和原始图像对血管提取的作用存在一定的差异，因此，为它们设置不同的权重因子后得到局部邻域的能量泛函为：
其中，参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子。
血管的提取问题在统计学上可以看成像素的分类问题，即将血管图像和原始图像作为两个相互独立的随机变量构建局部高斯拟合的能量泛函后，还需要对轮廓曲线的平滑性，曲线长度，以及局部轮廓细节进行约束，因此引入如下正则项：
其中，υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数(Heaviside function)，▽为梯度操作符，γ,κ为尺度参数。
进一步地，上述的方法中，使用变分水平集方法对能量泛函进行数学计算的步骤包括：根据变分法将能量泛函转化为偏微分方程的形式为:
其中，t表示时间参数，用来反映轮廓曲线随时间的变化情况；
δ(φ)为H(φ)函数针对自变量φ的导数；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
ω(x-y)为基于距离的高斯权重；
分别表示原始图像对应的局部灰度均值和标准差，分别表示增强的血管图像对应的局部灰度均值和标准差，它们均通过欧拉-拉格朗日方程确定；
参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数。
本发明的基于多特征高斯拟合的活动轮廓血管提取系统包括：
血管增强滤波模块，其对视网膜图像进行血管增强处理，即采用局部相位(localphase)血管增强算法将视网膜图像中血管对应的区域突显出来，降低复杂背景信息和像素灰度不均匀性对分割结果造成的干扰；其中，增强后的血管图像(vesselnessmap)可以作为一种新的特征图像用于血管提取，即在血管图像中，像素灰度值越高，其代表血管的可能性越大；反之则越小为图像背景的可能性越高大；
局部高斯拟合能量泛函模块，将血管图像和原始图像中的像素值作为两个相互独立的随机变量，构建一种基于二维高斯概率分布拟合的能量泛函活动轮廓模型，即在血管图像和原始视网膜图像中的像素值通过二维高斯分布统计模型模拟一定图像范围内的血管和背景区域内的图像统计特性，从而构建一个区域的活动轮廓模型。仅使用该能量泛函可能存在一定的不足，因此还需要通过一个正则约束项，惩罚血管轮廓曲线的不平滑性、不当的曲线长度、以及保持轮廓的局部细节；
变分求解框架模块，使用变分水平集方法对能量泛函进行数学计算，即血管轮廓以隐式水平集的形式表示，然后通过梯度下降流和欧拉-拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的血管轮廓。
进一步地，上述的系统中，血管增强滤波模块中，其采用的局部相位血管增强算法包括以下实施中的相关参数：设置中心频率为5π/7，带宽为2倍频(octave)，滤波范围为15×15，滤波方向为0°、45°、90°和135°，图像尺度系数为3，整合权重为3，正则系数为3。通过局部相位的血管增强可以显著突出血管所在的区域，同时减缓不相关背景信息的干扰。
进一步地，上述的系统中，局部高斯拟合能量泛函模块的能量泛函构建步骤如下：
在血管图像和原始视网膜图像中，对任一像素点局部邻域内的灰度概率分布通过二维高斯概率分布模型进行模拟：
P＝pi,x(I(y),V(y))＝pi,x(I(y))pi,x(V(y))；
其中，pi,x(I(y))为原始灰度图像中，中心像素点为x的邻域内，任一像素点y处的图像灰度I的概率分布，和为邻域内的灰度平均值和标准差；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
pi,x(V(y))，以及分别依次为增强后的血管图像V对应的概率分布、局部平均值和标准差；
血管的提取问题在统计学上可以看成像素的分类问题，即将血管所在区域的像素和背景区域的像素区分开，因此血管区域的像素点中，像素的灰度值应尽可能地接近于灰度平均值，即灰度标准差尽可能趋近于0；转化为概率模型后为：
Πypi,x(I(y))pi,x(V(y))→max；
将上式的最大化求解转化为最小化，并引入基于距离的高斯权重 其中，σ为高斯函数中的标准差参数，同时考虑到增强后的血管图像和原始图像对血管提取的作用存在一定的差异，因此，为它们设置不同的权重因子后得到局部邻域的能量泛函为：
其中，参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
将血管图像和原始图像作为两个相互独立的随机变量构建局部高斯拟合的能量泛函后，还需要对轮廓曲线的平滑性，曲线长度，以及局部轮廓细节进行约束，因此引入如下正则项：
其中，υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数。
进一步地，上述的系统中，使用变分水平集方法对能量泛函进行数学计算的步骤包括：根据变分法将能量泛函转化为偏微分方程的形式为:
其中，t表示时间参数，用来反映轮廓曲线随时间的变化情况；
δ(φ)为H(φ)函数针对自变量φ的导数；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
ω(x-y)为基于距离的高斯权重；
分别表示原始图像对应的局部灰度均值和标准差，分别表示增强的血管图像对应的局部灰度均值和标准差，它们均通过欧拉-拉格朗日方程确定；
参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数。
基于本发明的公开的视网膜图像数据(Digital Retinal Images for VesselExtraction,DRIVE)的仿真实验表明：本发明能够较为准确完整地提取出视网膜图像中各种不同粗细的血管组织，为眼部疾病的治疗提供有效的辅助，因此解决上述的技术问题。
附图说明
图1是本发明的基于多特征高斯拟合的活动轮廓血管提取方法的流程图。
图2是本发明的基于多特征高斯拟合的活动轮廓血管提取系统的结构框图。
图3A-D是本发明的仿真实验中的原始视网膜图像和血管增强后的血管图像。
图4A-B是本发明的仿真实验中提取的任选的一张图像结果图。
图5A-B是本发明的仿真实验中提取的不同正则项权重的分割结果图。
图6A-D是本发明的仿真实验中提取的不同活动轮廓算法的对比结果图。
具体实施方式
下面结合附图详细说明本发明的优选实施方式。
为了达到本发明的目的，如图1-2所示，在本发明的其中一些实施方式中，其提供了一种基于多特征高斯拟合的活动轮廓血管提取方法以及一种基于多特征高斯拟合的活动轮廓血管提取系统。
如图1-2所示，本发明的方法及系统基于多特征高斯拟合的活动轮廓血管提取算法，通过以下实施步骤实现：
步骤1，对输入视网膜图像进行血管增强处理：
步骤1a：对眼部视网膜图像进行多尺度多方向正交滤波，滤波器由相互正交且相位差为π/2的两个分滤波器组成(分别称为偶数滤波器和奇数滤波分量)，滤波结果以(a+bi)的复数形式构成一个复数图像数据，其实部和虚部分别对应于图像中的血管所在的区域以及血管的边界。具体求解公式为：e,o分别代表偶数和奇数滤波结果，且E,O分别为偶数和奇数滤波分量，m,j分别为滤波的尺度数和滤波方向，*代表卷积运算符，为复数单位；
步骤1b：对生成的复数图像数据进行多方向多尺度的整合运算并正则化处理，从而得到一个血管组织更加突出的血管增强结果。本文使用正则化处理复数图像的实部作为一种新的原图像特征(vesselness map)，用于血管提取算法中；
步骤2，局部高斯拟合的能量泛函：
将原图像中任选的一个像素点的灰度值和其对应的血管值作为两个相互独立的随机变量，然后在选定像素点的邻域内使用二维高斯概率分布模型量化这两个变量的统计特性。获得变量的分布特性后，以它们在目标区域和背景区域同时发生的概率的最大化为目标构建能量泛函；仅使用概率分布最大化构建的能量泛函是不足以获得满意结果的，因此，需要引入正则项用以惩罚轮廓曲线的不光滑性、不当的曲线长度、以及保持血管的局部细节信息。整合能量泛函和正则项后得到最终的活动轮廓模型，使用该模型进行血管提取，其获取的轮廓曲线将是短而平滑的，且能保持一定的局部细节；
步骤3，变分求解框架：
对最终生成的目标函数进行最大化通常转化为最小化的求解的过程，并使用变分法将其转化为偏微分方程的求解，即使用梯度下降流和欧拉-拉格朗日方程获得所需的一个关于时间的偏微分方程，并以迭代逼近的方式获取最终的血管轮廓。具体求解微分方程为：
步骤4，眼部血管轮廓提取：
迭代逼近的求解偏微分方程需要设置初始的轮廓曲线，并设定曲线内表示血管所在的区域，曲线外围背景区域。在微分方程的迭代计算下，轮廓曲线逐渐逼近所求血管的边界，从而将所需的血管提取出来。
下面介绍本发明的公开的视网膜图像数据的仿真实验：
(1)仿真条件：
本发明的仿真在Win7-64Intel(R)Core(TM)i3-2100CPU@3.10GHz 3.10GHz RAM6GB平台上的MATLAB 2013a软件上进行仿真模拟的，仿真数据选用荷兰乌UtrechtUniversity公开发表视网膜影像数据(Digital Retinal Images forVesselExtraction,DRIVE)；该数据分为训练数据和测试数据两部分，每部分有20幅图像，图像大小为768×584；此外，该套数据还提供两组人工提取的血管作为实际所需的血管组织；
(2)仿真内容与结果：
(2-1)仿真实验1：
本仿真实验使用数据集中的全部测试图像进行血管提取使用，验证算法的有效性，实验结果分别呈现在图3A-D和图4A-B中：从眼部视网膜图像中，可以明显看出复杂的背景信息，以及逐渐细化的血管走向，以及严重的灰度不均匀等导致血管与背景难以区分的问题。
图3A、B、C和D分别对应于原始视网膜图像区域，增强处理后的复数图像、复数实部图像和复数虚部图像；从图中可以看出，实部图像突出了血管所在区域，并且处在该区域的像素灰度值均大于零；而虚部则突出了血管的边界。
图4A-B分别表示本文算法结果以及其与手工分割间的性能对比图，从图4A中可以看出，算法能够较为准确完整地将不同粗细的血管提取出来，图4B则反映了算法分割结果与手工分割结果间不存在较为严重的差异，即分割算法能够到达较高的分割精度。
(2-2)仿真实验2：
本仿真实验通过设置不同的正则项权重系数，验证其对血管提取的影响程度，实验结果如图5A-B。
图5A-B分别为不同权重系数对应的分割结果，权重系数过大时，如图5A所示，分割结果将丢失掉大量的局部细节；权重系数过小时，如图5B所示，分割结果则将不相关的背景信息也提取出来(箭头所指出了算法提取的不相关的背景信息)。
图5A-B表明，正则项对血管提取的有效性有一定的影响，但其权重系数的设置要适当，否则将导致血管局部细节的丢失或者提取不相关的局部细节。
(2-2)仿真实验3：
在本仿真实验中，通过对比已有的基于区域的活动轮廓模型(即：CV model、LBFmodel和LGD model)与本文算法的分割结果，验证本发明的有效性。
使用测试数据进行的仿真实验中，实验结果通过敏感性(sensitivity,Se)、具体性(specificity,Sp)、准确性(accuracy,Acc)以及特征曲线面积(the area under areceiver operating characteristic curve,Auc)四个指标进行量化。它们分别表示如下：
其中，tp,fn,tn以及fp分别表示提取的真实血管像素位置(true positive)、错误提取的血管像素位置(false positive)、真实的背景像素位置(true negative)、以及错误的背景像素位置(false negative)。
Se和Sp均表示算法的有效性，前者为提取所需目标的有效性，后者则为提取背景信息的有效性。Acc反映了分割算法的一个综合性能。Auc是对Se和Sp的一个有效的折中。Se、Sp、Acc以及Auc的值越大，表示算法的分割性能越好。
仿真实验的对比结果参考图6A-D，其中，图6A、B、C和D分别对应于CV、LBF、LGDFmodels和本文方法的分割结果。从图中能够较为直接地观察到它们在分割性能上的差异，其中CV model的分割结果是四者中最为粗糙的，本文算法则是提取性能最优的。至于LBF，LGDF models和本文算法，它们间的差异可以通过图中圆环所在区域内的细小血管分割情况进行比较，即LBF不能提取指定区域内的细小血管，LGDF仅能够识别其中的一小部分，而本文算法能够提取其中的大部分。
对比结果表明，本文发明在血管提取性能上优于现存的几种活动轮廓模型。
对比实验的统计分析结果如下的表1所示，Se、Sp、Acc以及Auc四个评价指标的平均值和标准差数据参考表1。
从表1中可以看出：本文算法在四个评价指标上均能获得最大的平均值，其对应的标准差也相对由于其他三种分割算法。
因此，基于本发明的公开的视网膜图像数据(Digital Retinal ImagesforVessel Extraction,DRIVE)的仿真实验表明：本发明能够较为准确完整地提取出视网膜图像中各种不同粗细的血管组织，为眼部疾病的治疗提供有效的辅助，因此解决上述的技术问题。
以上所述的仅是本发明的优选实施方式，应当指出，对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。
Claims (6)
1.基于多特征高斯拟合的活动轮廓血管提取方法，其特征在于，包括：
步骤1，通过血管增强滤波模块对视网膜图像进行血管增强处理，即采用局部相位血管增强算法将视网膜图像中血管对应的区域突显出来，降低复杂背景信息和像素灰度不均匀性对分割结果造成的干扰；
步骤2，通过局部高斯拟合能量泛函模块将血管图像和原始图像中的像素值作为两个相互独立的随机变量，构建一种基于二维高斯概率分布拟合的能量泛函活动轮廓模型，其具体步骤如下：在血管图像和原始视网膜图像中，对任一像素点局部邻域内的灰度概率分布通过二维高斯概率分布模型进行模拟：
P＝pi,x(I(y),V(y))＝pi,x(I(y))pi,x(V(y))；
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>I</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mi>exp</mi>
<mo>(</mo>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>I</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
<mo>)</mo>
<mo>;</mo>
</mrow>
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>V</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mi>exp</mi>
<mo>(</mo>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>V</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
<mo>)</mo>
<mo>;</mo>
</mrow>
其中，pi,x(I(y))为原始灰度图像中，中心像素点为x的邻域内，任一像素点y处的图像灰度I的概率分布，和为邻域内的灰度平均值和标准差；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
pi，x(V(y))，以及分别依次为增强后的血管图像V对应的概率分布、局部平均值和标准差；
将血管所在区域的像素和背景区域的像素区分开，转化为概率模型后为：
Πypi,x(I(y))pi,x(V(y))→max；
将上式的最大化求解转化为最小化，并引入基于距离的高斯权重 其中，σ为高斯函数中的标准差参数，同时考虑到增强后的血管图像和原始图像对血管提取的作用存在一定的差异，因此，为它们设置不同的权重因子后得到局部邻域的能量泛函为：
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>E</mi>
<mi>x</mi>
</msub>
<mo>=</mo>
<mo>-</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mi>log</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>p</mi>
<mrow>
<mn>1</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>I</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
</msup>
<msub>
<mi>p</mi>
<mrow>
<mn>1</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>V</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
</msup>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>y</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mi>log</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>p</mi>
<mrow>
<mn>2</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>I</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
</msup>
<msub>
<mi>p</mi>
<mrow>
<mn>2</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>V</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
</msup>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>y</mi>
<mo>;</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中，参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
将血管图像和原始图像作为两个相互独立的随机变量构建局部高斯拟合的能量泛函后，还需要对轮廓曲线的平滑性，曲线长度，以及局部轮廓细节进行约束，因此引入如下正则项：
<mrow>
<mi>R</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>&upsi;</mi>
<mo>&Integral;</mo>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>|</mo>
<mi>d</mi>
<mi>x</mi>
<mo>+</mo>
<mi>&mu;</mi>
<mo>&Integral;</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<msup>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>&phi;</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>|</mo>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mi>d</mi>
<mi>x</mi>
<mo>+</mo>
<mi>&eta;</mi>
<mo>&Integral;</mo>
<mo>(</mo>
<mfrac>
<msup>
<mi>&gamma;</mi>
<mi>&kappa;</mi>
</msup>
<mrow>
<msup>
<mi>&gamma;</mi>
<mi>&kappa;</mi>
</msup>
<mo>+</mo>
<mi>&phi;</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mi>&kappa;</mi>
</msup>
</mrow>
</mfrac>
<mo>)</mo>
<mi>d</mi>
<mi>x</mi>
<mo>;</mo>
</mrow>
其中，υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数；
步骤3，通过变分求解框架模块，使用变分水平集方法对能量泛函进行数学计算，即血管轮廓以隐式水平集的形式表示，然后通过梯度下降流和欧拉-拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的血管轮廓。
2.根据权利要求1所述的基于多特征高斯拟合的活动轮廓血管提取方法，其特征在于，所述的局部相位血管增强算法包括以下实施中的相关参数：设置中心频率为5π/7，带宽为2倍频，滤波范围为15×15，滤波方向为0°、45°、90°和135°，图像尺度系数为3，整合权重为3，正则系数为3。
3.根据权利要求1所述的基于多特征高斯拟合的活动轮廓血管提取方法，其特征在于，所述步骤3中，使用变分水平集方法对能量泛函进行数学计算的具体步骤包括：根据变分法将能量泛函转化为偏微分方程的形式为:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>V</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>1</mn>
<mi>V</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>V</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>2</mn>
<mi>V</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>&upsi;</mi>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>i</mi>
<mi>v</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>&phi;</mi>
<mo>|</mo>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mo>&dtri;</mo>
<mn>2</mn>
</msup>
<mi>&phi;</mi>
<mo>-</mo>
<mi>d</mi>
<mi>i</mi>
<mi>v</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>&phi;</mi>
<mo>|</mo>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>&eta;</mi>
<mfrac>
<mrow>
<msup>
<mi>&kappa;&gamma;</mi>
<mi>&kappa;</mi>
</msup>
</mrow>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&gamma;</mi>
<mi>&kappa;</mi>
</msup>
<mo>+</mo>
<mi>&phi;</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mi>&kappa;</mi>
</msup>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mfrac>
<msup>
<mi>&phi;</mi>
<mrow>
<mo>(</mo>
<mi>&kappa;</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>;</mo>
</mrow>
<mrow>
<msubsup>
<mi>e</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<mi>I</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>y</mi>
<mo>;</mo>
</mrow>
<mrow>
<msubsup>
<mi>e</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<mi>V</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>y</mi>
<mo>;</mo>
</mrow>
其中，t表示时间参数，用来反映轮廓曲线随时间的变化情况；
δ(φ)为H(φ)函数针对自变量φ的导数；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
ω(x-y)为基于距离的高斯权重；
分别表示原始图像对应的局部灰度均值和标准差，分别表示增强的血管图像对应的局部灰度均值和标准差，它们均通过欧拉-拉格朗日方程确定；
参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数。
4.基于多特征高斯拟合的活动轮廓血管提取系统，其特征在于，包括：
血管增强滤波模块，其对视网膜图像进行血管增强处理，即采用局部相位血管增强算法将视网膜图像中血管对应的区域突显出来，降低复杂背景信息和像素灰度不均匀性对分割结果造成的干扰；
局部高斯拟合能量泛函模块，将血管图像和原始图像中的像素值作为两个相互独立的随机变量，构建一种基于二维高斯概率分布拟合的能量泛函活动轮廓模型，其具体步骤如下：在血管图像和原始视网膜图像中，对任一像素点局部邻域内的灰度概率分布通过二维高斯概率分布模型进行模拟：
P＝pi,x(I(y),V(y))＝pi,x(I(y))pi,x(V(y))；
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>I</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mi>exp</mi>
<mo>(</mo>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<mi>I</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
<mo>)</mo>
<mo>;</mo>
</mrow>
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>V</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mi>exp</mi>
<mo>(</mo>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<mi>V</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
<mo>)</mo>
<mo>;</mo>
</mrow>
其中，pi,x(I(y))为原始灰度图像中，中心像素点为x的邻域内，任一像素点y处的图像灰度I的概率分布，和为邻域内的灰度平均值和标准差；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
pi，x(V(y))，以及分别依次为增强后的血管图像V对应的概率分布、局部平均值和标准差；
将血管所在区域的像素和背景区域的像素区分开，转化为概率模型后为：
∏ypi,x(I(y))pi,x(V(y))→max；
将上式的最大化求解转化为最小化，并引入基于距离的高斯权重 其中，σ为高斯函数中的标准差参数，同时考虑到增强后的血管图像和原始图像对血管提取的作用存在一定的差异，因此，为它们设置不同的权重因子后得到局部邻域的能量泛函为：
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>E</mi>
<mi>x</mi>
</msub>
<mo>=</mo>
<mo>-</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mi>log</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>p</mi>
<mrow>
<mn>1</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>I</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
</msup>
<msub>
<mi>p</mi>
<mrow>
<mn>1</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>V</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
</msup>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>y</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mi>log</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>p</mi>
<mrow>
<mn>2</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>I</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
</msup>
<msub>
<mi>p</mi>
<mrow>
<mn>2</mn>
<mo>,</mo>
<mi>x</mi>
</mrow>
</msub>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>V</mi>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
</msup>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>y</mi>
<mo>;</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中，参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
将血管图像和原始图像作为两个相互独立的随机变量构建局部高斯拟合的能量泛函后，还需要对轮廓曲线的平滑性，曲线长度，以及局部轮廓细节进行约束，因此引入如下正则项：
<mrow>
<mi>R</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>&upsi;</mi>
<mo>&Integral;</mo>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>H</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>|</mo>
<mi>d</mi>
<mi>x</mi>
<mo>+</mo>
<mi>&mu;</mi>
<mo>&Integral;</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<msup>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>&phi;</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>|</mo>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mi>d</mi>
<mi>x</mi>
<mo>+</mo>
<mi>&eta;</mi>
<mo>&Integral;</mo>
<mo>(</mo>
<mfrac>
<msup>
<mi>&gamma;</mi>
<mi>&kappa;</mi>
</msup>
<mrow>
<msup>
<mi>&gamma;</mi>
<mi>&kappa;</mi>
</msup>
<mo>+</mo>
<mi>&phi;</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mi>&kappa;</mi>
</msup>
</mrow>
</mfrac>
<mo>)</mo>
<mi>d</mi>
<mi>x</mi>
<mo>;</mo>
</mrow>
其中，υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数；
变分求解框架模块，使用变分水平集方法对能量泛函进行数学计算，即血管轮廓以隐式水平集的形式表示，然后通过梯度下降流和欧拉-拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的血管轮廓。
5.根据权利要求4所述的基于多特征高斯拟合的活动轮廓血管提取系统，其特征在于，所述的血管增强滤波模块中，其采用的局部相位血管增强算法包括以下实施中的相关参数：设置中心频率为5π/7，带宽为2倍频，滤波范围为15×15，滤波方向为0°、45°、90°和135°，图像尺度系数为3，整合权重为3，正则系数为3。
6.根据权利要求4所述的基于多特征高斯拟合的活动轮廓血管提取系统，其特征在于，所述的变分求解框架模块中使用变分水平集方法对能量泛函进行数学计算的步骤包括：根据变分法将能量泛函转化为偏微分方程的形式为:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>1</mn>
<mi>I</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>2</mn>
<mi>I</mi>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>1</mn>
<mi>V</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>1</mn>
<mi>V</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&lambda;</mi>
<mn>2</mn>
<mi>V</mi>
</msubsup>
<msubsup>
<mi>e</mi>
<mn>2</mn>
<mi>V</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>&upsi;</mi>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<mi>&phi;</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>i</mi>
<mi>v</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>&phi;</mi>
<mo>|</mo>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mo>&dtri;</mo>
<mn>2</mn>
</msup>
<mi>&phi;</mi>
<mo>-</mo>
<mi>d</mi>
<mi>i</mi>
<mi>v</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mo>&dtri;</mo>
<mi>&phi;</mi>
</mrow>
<mrow>
<mo>|</mo>
<mo>&dtri;</mo>
<mi>&phi;</mi>
<mo>|</mo>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>&eta;</mi>
<mfrac>
<mrow>
<msup>
<mi>&kappa;&gamma;</mi>
<mi>&kappa;</mi>
</msup>
</mrow>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&gamma;</mi>
<mi>&kappa;</mi>
</msup>
<mo>+</mo>
<mi>&phi;</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mi>&kappa;</mi>
</msup>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mfrac>
<msup>
<mi>&phi;</mi>
<mrow>
<mo>(</mo>
<mi>&kappa;</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>;</mo>
</mrow>
<mrow>
<msubsup>
<mi>e</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<mi>I</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>I</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>y</mi>
<mo>;</mo>
</mrow>
<mrow>
<msubsup>
<mi>e</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>&Integral;</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msqrt>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</msqrt>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>u</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<mi>V</mi>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>i</mi>
<mi>V</mi>
</msubsup>
<msup>
<mrow>
<mo>(</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>y</mi>
<mo>;</mo>
</mrow>
其中，t表示时间参数，用来反映轮廓曲线随时间的变化情况；
δ(φ)为H(φ)函数针对自变量φ的导数；
i表示轮廓曲线的内外，当i＝1表示像素点位于曲线内，当i＝2时表示像素点位于曲线外；
ω(x-y)为基于距离的高斯权重；
分别表示原始图像对应的局部灰度均值和标准差，分别表示增强的血管图像对应的局部灰度均值和标准差，它们均通过欧拉-拉格朗日方程确定；
参数和为轮廓曲线内原始图像和增强血管图像对应的权重因子；参数和为轮廓曲线外原始图像和增强血管图像对应的权重因子；
υ,μ,η是三种不同的权重系数，分别控制着轮廓曲线的平滑性，曲线长度以及曲线的局部拓扑；φ(x)为零水平集函数，H(φ(x))为以φ(x)为自变量的单位阶跃函数，▽为梯度操作符，γ,κ为尺度参数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510504854.8A CN105160660B (zh) | 2015-08-17 | 2015-08-17 | 基于多特征高斯拟合的活动轮廓血管提取方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510504854.8A CN105160660B (zh) | 2015-08-17 | 2015-08-17 | 基于多特征高斯拟合的活动轮廓血管提取方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105160660A CN105160660A (zh) | 2015-12-16 |
| CN105160660B true CN105160660B (zh) | 2017-12-01 |
Family
ID=54801502
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510504854.8A Active CN105160660B (zh) | 2015-08-17 | 2015-08-17 | 基于多特征高斯拟合的活动轮廓血管提取方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105160660B (zh) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106570882B (zh) * | 2016-10-28 | 2019-08-20 | 中国科学院苏州生物医学工程技术研究所 | 混合高斯分布模型的活动轮廓图像分割方法 |
| CN106548478B (zh) * | 2016-10-28 | 2019-06-28 | 中国科学院苏州生物医学工程技术研究所 | 基于局部拟合图像的活动轮廓图像分割方法 |
| CN106530314A (zh) * | 2016-12-21 | 2017-03-22 | 中国科学院合肥物质科学研究院 | 一种多尺度局部统计主动轮廓模型(lsacm)水平集图像分割方法 |
| CN106934815A (zh) * | 2017-02-27 | 2017-07-07 | 南京理工大学 | 基于混合区域的活动轮廓模型图像分割方法 |
| CN108550180B (zh) * | 2018-03-09 | 2021-11-19 | 南京信息工程大学 | 基于内点集域约束及高斯过程参数优化的血管建模方法 |
| CN108460782A (zh) * | 2018-04-03 | 2018-08-28 | 信阳农林学院 | 基于主动轮廓模型的生猪耳部图像分割方法 |
| CN110223771B (zh) * | 2019-05-29 | 2023-05-30 | 安徽医科大学第一附属医院 | 基于nhpp的消化内科电子数据分析方法 |
| CN110648340B (zh) * | 2019-09-29 | 2023-03-17 | 惠州学院 | 一种基于二进制及水平集处理图像的方法及装置 |
| CN110751634A (zh) * | 2019-10-11 | 2020-02-04 | 北京致远慧图科技有限公司 | 一种视杯视盘分割模型的确定方法、装置及存储介质 |
| CN110706225B (zh) * | 2019-10-14 | 2020-09-04 | 山东省肿瘤防治研究院（山东省肿瘤医院） | 基于人工智能的肿瘤识别系统 |
| CN115187598B (zh) * | 2022-09-09 | 2023-02-03 | 天津远景科技服务有限公司 | 血管造影图像的处理方法、装置、系统、设备及介质 |
| CN115457054A (zh) * | 2022-09-13 | 2022-12-09 | 深圳先进技术研究院 | 图像分割方法、装置、设备及可读存储介质 |
| CN116740768B (zh) * | 2023-08-11 | 2023-10-20 | 南京诺源医疗器械有限公司 | 基于鼻颅镜的导航可视化方法、系统、设备及存储介质 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102982547A (zh) * | 2012-11-29 | 2013-03-20 | 北京师范大学 | 自动初始化的局域活动轮廓模型心脑血管分割方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8755576B2 (en) * | 2011-09-09 | 2014-06-17 | Calgary Scientific Inc. | Determining contours of a vessel using an active contouring model |
-
2015
- 2015-08-17 CN CN201510504854.8A patent/CN105160660B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102982547A (zh) * | 2012-11-29 | 2013-03-20 | 北京师范大学 | 自动初始化的局域活动轮廓模型心脑血管分割方法 |
Non-Patent Citations (4)
| Title |
|---|
| Active contours driven by local Gaussian distribution fitting energy;Li Wang et al.;《Signal Processing》;20091231;第89卷(第12期);第2435-2447页 * |
| Minimization of Region-Scalable Fitting Energy for Image Segmentation;Chunming Li et al.;《IEEE Transactions on Image Processing》;20081031;第17卷(第10期);第1940-1949页 * |
| Numerical Conditioning Problems and Solutions for Nonparametric i.i.d. Statistical Active Contours;Hao Wu et al.;《IEEE Transactions on Pattern Analysis and Machine Intelligence》;20130630;第35卷(第6期);第1298-1311页 * |
| 改进C-V分割算法在多光谱成像仪中的应用;张艳超 等;《中国光学》;20150228;第8卷(第1期);第68-73页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105160660A (zh) | 2015-12-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105160660B (zh) | 基于多特征高斯拟合的活动轮廓血管提取方法及系统 | |
| Aguirre-Ramos et al. | Blood vessel segmentation in retinal fundus images using Gabor filters, fractional derivatives, and Expectation Maximization | |
| CN109166124B (zh) | 一种基于连通区域的视网膜血管形态量化方法 | |
| CN107610087B (zh) | 一种基于深度学习的舌苔自动分割方法 | |
| CN104809740B (zh) | 基于svm与弹性区域生长的膝软骨图像自动分割方法 | |
| CN108257135A (zh) | 基于深度学习方法解读医学图像特征的辅助诊断系统 | |
| CN104899876B (zh) | 一种基于自适应高斯差分的眼底图像血管分割方法 | |
| CN109345538A (zh) | 一种基于卷积神经网络的视网膜血管分割方法 | |
| CN107248161A (zh) | 一种多特征融合的有监督视网膜血管提取方法 | |
| CN106934816A (zh) | 一种基于elm的眼底图像视网膜血管分割方法 | |
| CN110276763B (zh) | 一种基于可信度和深度学习的视网膜血管分割图生成方法 | |
| CN106408017A (zh) | 基于深度学习的超生颈动脉内中膜厚度测量装置和方法 | |
| CN104809480A (zh) | 一种基于分类回归树和AdaBoost的眼底图像视网膜血管分割方法 | |
| CN106485721A (zh) | 从光学相干断层图像获取视网膜结构的方法及其系统 | |
| Zhou et al. | An iterative speckle filtering algorithm for ultrasound images based on bayesian nonlocal means filter model | |
| CN108053398A (zh) | 一种半监督特征学习的黑色素瘤自动检测方法 | |
| CN103942780B (zh) | 基于改进模糊连接度算法的丘脑及其子结构分割方法 | |
| Samiappan et al. | Classification of carotid artery abnormalities in ultrasound images using an artificial neural classifier. | |
| CN101667289A (zh) | 基于nsct特征提取和监督分类的视网膜图像分割方法 | |
| CN111862009A (zh) | 一种眼底oct图像的分类方法及计算机可读存储介质 | |
| Chen et al. | Retina image vessel segmentation using a hybrid CGLI level set method | |
| Yang et al. | Classification of diabetic retinopathy severity based on GCA attention mechanism | |
| CN106934815A (zh) | 基于混合区域的活动轮廓模型图像分割方法 | |
| Parvathi et al. | Automatic drusen detection from colour retinal images | |
| Lyu et al. | Deep tessellated retinal image detection using Convolutional Neural Networks |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |