CN111028241B - 一种多尺度血管增强的水平集分割系统与方法 - Google Patents

Abstract

本发明提供一种多尺度血管增强的水平集分割系统与方法，涉及图像处理技术领域。本发明包括输入模块、预处理模块、血管增强模块、水平集模块以及输出模块；本发明针对以往血管图像中，待分割的血管与背景区域往往呈现出灰度差异小，灰度分布不均匀的现象，克服了血管图像中血管与背景灰度差异不明显对分割精度的影响，在多尺度的血管增强下，增加了水平集分割方法对不同粗细血管的分割能力，实现微细血管的精准分割。

Description

一种多尺度血管增强的水平集分割系统与方法

技术领域

本发明涉及图像处理技术领域，尤其涉及一种多尺度血管增强的水平集分割系统与方法。

背景技术

血管分割是现代医学影像辅助诊断系统中的关键步骤之一，通过血管分割可以更直观地研究血管特征，进而对相关疾病的预防，诊断和治疗提出重要参考依据。水平集方法是广泛应用于血管分割的一类方法，主要思想是用平滑封闭的曲线标记感兴趣边界，通过不断迭代来实现分割。但是传统的水平集方法很大程度上依赖图像中背景和血管之间的灰度差异，而在诸多场景中，待分割的血管与背景区域往往呈现出灰度差异小，灰度分布不均匀的现象。本发明提出了一种多尺度血管增强的水平集分割系统与方法，解决血管粗细粒度以及图像质量给血管分割精度带来的挑战。

发明内容

针对现有技术存在的问题，本发明提供一种多尺度血管增强的水平集分割系统与方法，

本发明所采取的技术方案是：一种多尺度血管增强的水平集分割系统，包括输入模块、预处理模块、血管增强模块、水平集模块以及输出模块；

所述输入模块，用于接收用户输入的眼底血管图像，然后将其传递到预处理模块；

所述预处理模块，用于对接收的眼底血管图像进行扩展，扩展的结果传递到血管增强模块；

所述血管增强模块，用于接收预处理后的图片，生成图像M，并将处理后的图片M传递到水平集模块

所述水平集模块用于对图像M采用水平集算法进行血管分割，并将分割结果传递到输出模块；

所述输出模块用于对分割结果图像进行输出显示。

另一方面，本发明提供多尺度血管增强的水平集分割方法，通过前述多尺度血管增强的水平集分割系统实现，包括以下步骤：

步骤1：输入模块接收待分割血管图像I，将图像I传递到预处理模块；

步骤2：在预处理模块中，将图像I的定义域定义为二维实数空间子集Ω；

步骤2.1：若I是彩色图像，则令J＝αI(R)+βI(G)+γI(B)；其中I(R),I(G),I(B)分别为图像I的三通道，其中J是通过改变权重后得到的新图像，α，β，γ分别是各通道的权重；

若I不是彩色图像，则J＝I；

步骤2.2：构造无边缘背景的眼底图像M，其分辨率与图像I的分辨率相同，若待分割区域充满Ω，则M＝J，若待分割区域没有充满Ω，对图像内的有效边界区域进行镜像扩展使整个图像内的黑色填充区被覆盖，扩展后的图像为P，使M＝P；

所述镜像扩展方法为：

步骤2.2.1：确定阈值Q将J划分为前景R和背景B；

步骤2.2.2：将R进行腐蚀操作得到R₁，构造前景区域边界图像R₂，使R₂＝R-R₁；

步骤2.2.3：设点W是背景B上一个点，遍历R₂中的每一个点，使得W与R_p(R_p∈R₂)的距离最小，即min(dist(W,R_p))，以R_p为中点与W做中心对称，可得到R中一个对应像素点，将该点的值赋给当前W；

步骤2.2.4：对背景B中的每一个点重复步骤2.2.3操作，直到每一个点被赋值，即完成图像扩展得到图像M；

步骤2.2.5：将图像M传递到血管增强模块。

步骤3：在血管增强模块中，设(u，v)为图像M中任意一点，用多尺度Hessian矩阵构造Frangi滤波器获得该点增强值F(u,v)，遍历图像M获取每个点的增强值即得到血管增强图像F；

步骤3.1：设定参数β，系数c，系数σ_g，尺度步长△σ，其中σ_g的尺度范围为[σ_min，σ_max]，并进行初始化令σ_g＝σ_min，且σ_min＜σ_max，F(u，v)＝0；

步骤3.2：构造高斯-拉普拉斯算子，即LOG算子，其中，基于尺度σ_g的高斯函数为：

可得出LOG算子如下：

其中，G是高斯函数，σg是高斯函数的尺度，x，y为坐标点横纵坐标，Gxx，Gyy，Gxy分别是高斯函数的二阶偏导数

步骤3.3：将图像M与LOG算子卷积构造Hessian矩阵，即矩阵H的公式为

其中，

r_xx＝G_xx*M，r_xy＝G_xy*M，r_yy＝G_yy*M，

其中*为卷积符号

步骤3.4：计算步骤3.3中得到的Hessian矩阵的两个特征值λ₁，λ₂

其中，K＝(r_xx+r_yy)/2，Q＝(r_xxr_yy-r_xyr_yx)/2，管状血管的Hessian矩阵的特征值满足|λ₁|≈0，|λ₁|＜＜|λ₂|

步骤3.5：由步骤3.4求得的Hessian矩阵特征值得到(u,v)的增强值F(u,v)：

其中，β和c都是公式参数，都是常数系数；

若F(u,v)＜f，则F(u,v)＝f，判断σ_g是否达到设定最大尺度σ_max，若未达到，则更新σ_g＝σ_g+△σ，转到步骤3.2，否则迭代结束，得到增强值F(u,v)；

步骤3.6：将管增强图像F传递到水平集模块中；

步骤4：在水平集模块中，对图像F进行血管分割，具体步骤如下：

步骤4.1：对用于划分图像区域的水平集函数进行初始化，在待分割图像F上对零水平集分割曲线进行初始化设定，根据设定的零水平集生成高一维的水平集函数φ；

步骤4.2：建立图像水平集分割的能量泛函，并设定水平集方法中的参数λ₁，λ₂，水平集分割曲线平滑项系数μ，面积平滑项系数v以及水平集函数正则项系数γ；

所述能量泛函如下：

E＝E_en(c₁,c₂,φ)+μLength(φ)+νArea(φ)+γP(φ)，

其中，E_en(c₁,c₂,φ)＝λ₁∫_Ω|F-c₁|²H(φ(u,v))dudv+λ₂∫_Ω|F-c₂|²(1-H(φ(u,v)))dudv为图像数据项，和Area(φ)＝∫H(φ(u,v))dudv分别是分割对象的长度约束项和面积约束项，/>是正则化项，H(φ)是φ的Heaviside函数；

步骤4.3：设定最大迭代次数J和图像水平集分割的能量泛函E的收敛阈值ε；

步骤4.4：更新图像水平集函数φ^j、以及H(φ^j)，1-H(φ^j)，得到第j次迭代的能量泛函E^j，若当前迭代次数j未达到最大迭代次数J或者能量泛函E相邻两次的迭代结果差值大于提前设定的收敛阈值ε，则迭代次数加1，继续更新；否则达到迭代终止条件，停止迭代，得到最后的分割结果；

其中，

所述更新图像水平集函数φ^j的公式为：

其中，△t为时间步长，为迭代更新的水平集函数的时间偏导数，所述时间偏导数公式为：

其中，div表示散度，表示梯度；

步骤4.5：将步骤4.4中得到的分割结果作用步骤3.1中的感兴趣区域R，得到最后的分割结果；

步骤5：输出模块得到最后结果。

采用上述技术方案所产生的有益效果在于：

本发明提出了一种多尺度血管增强后的水平集分割系统及方法，克服了血管图像中血管与背景灰度差异不明显对分割精度的影响，在多尺度的血管增强下，增加了水平集分割方法对不同粗细血管的分割能力，实现微细血管的精准分割。

附图说明

图1为多尺度血管增强后的水平集分割系统模块图；

图2为实施例中多尺度血管加强后的水平集分割方法流程图；

图3为分割区域没有充满整个图像区域的示例图；

图4为实施例中的彩色眼底血管图像原图；

其中，图(a)为实施例1中的彩色眼底血管图像原图；

图(b)为实施例中改变权重后的灰度扩展图像；

图(c)为实施例中的血管增强后的眼底图；

图(d)为实施例中最后的分割结果；

图5为实施例中设定的初始零水平集位置，及其不同迭代次数下的分割结果图；

其中，图(a)为实施例中设定的初始水平集位置；

图(b)为实施例中迭代10次后的水平集分割结果；

图(c)为实施例中迭代30次后的水平集分割结果；

图(d)为实施例中迭代50次后的水平集分割结果。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

一种多尺度血管增强的水平集分割系统，如图1所述，包括输入模块、预处理模块、血管增强模块、水平集模块以及输出模块；

所述输入模块，用于接收用户输入的眼底血管图像，然后将其传递到预处理模块；

所述预处理模块，用于对接收的眼底血管图像进行扩展，扩展的结果传递到血管增强模块；

所述血管增强模块，用于接收预处理后的图片，生成图像M，并将处理后的图片M传递到水平集模块

所述水平集模块用于对图像M采用水平集算法进行血管分割，并将分割结果传递到输出模块；

所述输出模块用于对分割结果图像进行输出显示。

另一方面，本发明提供多尺度血管增强的水平集分割方法，通过前述多尺度血管增强的水平集分割系统实现，如图2所示，包括以下步骤：

步骤1：输入模块接收待分割血管图像I，如图4(a)所示，将图像I传递到预处理模块；

步骤2：在预处理模块中，将图像I的定义域定义为二维实数空间子集Ω；

步骤2.1：若I是彩色图像，则令J＝αI(R)+βI(G)+γI(B)；其中I(R),I(G),I(B)分别为图像I的三通道，其中J是通过改变权重后得到的新图像，α，β，γ分别是各通道的权重；

若I不是彩色图像，则J＝I，如图4(b)所示；

步骤2.2：构造无边缘背景的眼底图像M，其分辨率与图像I的分辨率相同，若待分割区域充满Ω，则M＝J，若待分割区域没有充满Ω，对图像内的有效边界区域进行镜像扩展使整个图像内的黑色填充区被覆盖，如图3所示，扩展后的图像为P，使M＝P；

所述镜像扩展方法为：

步骤2.2.1：确定阈值Q将J划分为前景R和背景B；

步骤2.2.2：将R进行腐蚀操作得到R₁，构造前景区域边界图像R₂，使R₂＝R-R₁；

步骤2.2.3：设点W是背景B上一个点，遍历R₂中的每一个点，使得W与R_p(R_p∈R₂)的距离最小，即min(dist(W,R_p))，以R_p为中点与W做中心对称，可得到R中一个对应像素点，将该点的值赋给当前W；

步骤2.2.4：对背景B中的每一个点重复步骤2.2.3操作，直到每一个点被赋值，即完成图像扩展得到图像M；

步骤2.2.5：将图像M传递到血管增强模块。

步骤3：在血管增强模块中，设(u，v)为图像M中任意一点，用多尺度Hessian矩阵构造Frangi滤波器获得该点增强值F(u,v)，遍历图像M获取每个点的增强值即得到血管增强图像F，如图4(c)所示；

步骤3.1：设定参数β，系数c，系数σ_g，尺度步长△σ，其中σ_g的尺度范围为[σ_min，σ_max]，并进行初始化令σ_g＝σ_min，且σ_min＜σ_max，F(u，v)＝0；

步骤3.2：构造高斯-拉普拉斯算子，即LOG算子，其中，基于尺度σ_g的高斯函数为：

可得出LOG算子如下：

其中，G是高斯函数，σg是高斯函数的尺度，x，y为坐标点横纵坐标，Gxx，Gyy，Gxy分别是高斯函数的二阶偏导数；

步骤3.3：将图像M与LOG算子卷积构造Hessian矩阵，即矩阵H的公式为

其中，

r_xx＝G_xx*M，r_xy＝G_xy*M，r_yy＝G_yy*M，

其中*为卷积符号

步骤3.4：计算步骤3.3中得到的Hessian矩阵的两个特征值λ₁，λ₂

其中，K＝(r_xx+r_yy)/2，Q＝(r_xxr_yy-r_xyr_yx)/2，管状血管的Hessian矩阵的特征值满足|λ₁|≈0，|λ₁|＜＜|λ₂|

步骤3.5：由步骤3.4求得的Hessian矩阵特征值得到(u,v)的增强值F(u,v)：

其中，β和c都是公式参数，都是常数系数；

若F(u,v)＜f，则F(u,v)＝f，判断σ_g是否达到设定最大尺度σ_max，若未达到，则更新σ_g＝σ_g+△σ，转到步骤3.2，否则迭代结束，得到增强值F(u,v)；

步骤3.6：将管增强图像F传递到水平集模块中；

步骤4：在水平集模块中，对图像F进行血管分割，具体步骤如下：

步骤4.1：对用于划分图像区域的水平集函数进行初始化，在待分割图像F上对零水平集分割曲线进行初始化设定，根据设定的零水平集生成高一维的水平集函数φ；

步骤4.2：建立图像水平集分割的能量泛函，并设定水平集方法中的参数λ₁，λ₂，水平集分割曲线平滑项系数μ，面积平滑项系数v以及水平集函数正则项系数γ；

所述能量泛函如下：

E＝E_en(c₁,c₂,φ)+μLength(φ)+νArea(φ)+γP(φ)，

其中，E_en(c₁,c₂,φ)＝λ₁∫_Ω|F-c₁|²H(φ(u,v))dudv+λ₂∫_Ω|F-c₂|²(1-H(φ(u,v)))dudv为图像数据项，和Area(φ)＝∫H(φ(u,v))dudv分别是分割对象的长度约束项和面积约束项，/>是正则化项，H(φ)是φ的Heaviside函数；

步骤4.3：设定最大迭代次数J和图像水平集分割的能量泛函E的收敛阈值ε；

步骤4.4：更新图像水平集函数φ^j、以及H(φ^j)，1-H(φ^j)，得到第j次迭代的能量泛函E^j，若当前迭代次数j未达到最大迭代次数J或者能量泛函E相邻两次的迭代结果差值大于提前设定的收敛阈值ε，则迭代次数加1，继续更新；否则达到迭代终止条件，停止迭代，得到最后的分割结果，其迭代过程如图5中图5(a)-(d)所示；

其中，

所述更新图像水平集函数φ^j的公式为：

其中，△t为时间步长，为迭代更新的水平集函数的时间偏导数，所述时间偏导数公式为：

其中，div表示散度，表示梯度；

步骤4.5：将步骤4.4中得到的分割结果作用步骤3.1中的感兴趣区域R，得到最后的分割结果；

步骤5：输出模块得到最后结果，如图4(d)所示。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (1)

1.一种多尺度血管增强的水平集分割系统，其特征在于：包括输入模块、预处理模块、血管增强模块、水平集模块以及输出模块；

所述输入模块，用于接收用户输入的血管图像，然后将其传递到预处理模块；

所述预处理模块，用于对接收的血管图像进行扩展，扩展的结果传递到血管增强模块；

所述血管增强模块，用于接收预处理后的图片，生成图像M，并将处理后的图片M传递到水平集模块；

所述水平集模块用于对图像M采用水平集算法进行血管分割，并将分割结果传递到输出模块；

所述输出模块用于对分割结果图像进行输出显示；

所述一种多尺度血管增强的水平集分割系统，用于实现多尺度血管增强的水平集分割方法，包括以下步骤：

步骤1：输入模块接收待分割血管图像I，将图像I传递到预处理模块；

步骤2：在预处理模块中，将图像I的定义域定义为二维实数空间子集Ω；

步骤2.1：若I是彩色图像，则令J＝αI(R)+βI(G)+γI(B)；其中I(R),I(G),I(B)分别为图像I的三通道，其中J是通过改变权重后得到的新图像，α，β，γ分别是各通道的权重；

若I不是彩色图像，则J＝I；

步骤2.2：构造无边缘背景的眼底图像M，其分辨率与图像I的分辨率相同，若待分割区域充满Ω，则M＝J，若待分割区域没有充满Ω，对图像内的有效边界区域进行镜像扩展使整个图像内的黑色填充区被覆盖，扩展后的图像为P，使M＝P；

所述镜像扩展方法为：

步骤2.2.1：确定阈值Q将J划分为前景R和背景B；

步骤2.2.2：将R进行腐蚀操作得到R₁，构造前景区域边界图像R₂，使R₂＝R-R₁；

步骤2.2.3：设点W是背景B上一个点，遍历R₂中的每一个点，使得W与R_p的距离最小，R_p∈R₂，即min(dist(W,R_p))，以R_p为中点与W做中心对称，可得到R中一个对应像素点，将该点的值赋给当前W；

步骤2.2.4：对背景B中的每一个点重复步骤2.2.3操作，直到每一个点被赋值，即完成图像扩展得到图像M；

步骤2.2.5：将图像M传递到血管增强模块；

步骤3：在血管增强模块中，设(u，v)为图像M中任意一点，用多尺度Hessian矩阵构造Frangi滤波器获得该点增强值F(u,v)，遍历图像M获取每个点的增强值即得到血管增强图像F；

步骤3中所述用多尺度Hessian矩阵构造Frangi滤波器获得该点增强值F(u,v)的具体步骤为：

步骤3.1：设定参数β，系数c，系数σ_g，尺度步长△σ，其中σ_g的尺度范围为[σ_min，σ_max]，并进行初始化令σ_g＝σ_min，且σ_min＜σ_max，F(u，v)＝0；

步骤3.2：构造高斯-拉普拉斯算子，即LOG算子，其中，基于尺度σ_g的高斯函数为：

可得出LOG算子如下：

其中，G是高斯函数，σ_g是高斯函数的尺度，x，y为坐标点横纵坐标，G_xx，G_yy，G_xy分别是高斯函数的二阶偏导数；

步骤3.3：将图像M与LOG算子卷积构造Hessian矩阵，即矩阵H的公式为

其中，

r_xx＝G_xx*M，r_xy＝G_xy*M，r_yy＝G_yy*M，

其中*为卷积符号

步骤3.4：计算步骤3.3中得到的Hessian矩阵的两个特征值λ₁，λ₂

其中，K＝(r_xx+r_yy)/2，Q＝(r_xxr_yy-r_xyr_yx)/2，管状血管的Hessian矩阵的特征值满足|λ₁|≈0，|λ₁|＜＜|λ₂|

步骤3.5：由步骤3.4求得的Hessian矩阵特征值得到(u,v)的增强值F(u,v)：

其中，β和c都是公式参数，都是常数系数；

若F(u,v)＜f，则F(u,v)＝f，判断σ_g是否达到设定最大尺度σ_max，若未达到，则更新σ_g＝σ_g+△σ，转到步骤3.2，否则迭代结束，得到增强值F(u,v)；

步骤3.6：将血管增强图像F传递到水平集模块中；

步骤4：在水平集模块中，对图像F进行血管分割，具体步骤如下：

步骤4.1：对用于划分图像区域的水平集函数进行初始化，在待分割图像F上对零水平集分割曲线进行初始化设定，根据设定的零水平集生成高一维的水平集函数φ；

步骤4.2：建立图像水平集分割的能量泛函，并设定水平集方法中的参数λ₁，λ₂，水平集分割曲线平滑项系数μ，面积平滑项系数v以及水平集函数正则项系数γ；

所述能量泛函如下：

E＝E_en(c₁,c₂,φ)+μLength(φ)+νArea(φ)+γP(φ)，

其中，E_en(c₁,c₂,φ)＝λ₁∫_Ω|F-c₁|²H(φ(u,v))dudv+λ₂∫_Ω|F-c₂|²(1-H(φ(u,v)))dudv为图像数据项，Length(φ)＝∫|▽H(φ(u,v))|dudv和Area(φ)＝∫H(φ(u,v))dudv分别是分割对象的长度约束项和面积约束项，是正则化项，H(φ)是φ的Heaviside函数；

步骤4.3：设定最大迭代次数J和图像水平集分割的能量泛函E的收敛阈值ε；

步骤4.4：更新图像水平集函数φ^j、以及H(φ^j)，1-H(φ^j)，得到第j次迭代的能量泛函E^j，若当前迭代次数j未达到最大迭代次数J或者能量泛函E相邻两次的迭代结果差值大于提前设定的收敛阈值ε，则迭代次数加1，继续更新；否则达到迭代终止条件，停止迭代，得到最后的分割结果；

其中，

所述更新图像水平集函数φ^j的公式为：

其中，△t为时间步长，为迭代更新的水平集函数的时间偏导数，所述时间偏导数公式为：

其中，div表示散度，▽表示梯度；

步骤4.5：将步骤4.4中得到的分割结果作用步骤3.1中的感兴趣区域R，得到最后的分割结果；

步骤5：输出模块得到最后结果。