一种基于统计学模型和多尺度滤波的脑血管图像分割方法技术

本发明专利技术提出了一种基于统计学模型和多尺度滤波的脑血管图像分割方法，属于医学计算及图像处理技术领域。所述分割方法将统计模型和多尺度滤波的相结合，以在高层与低层模型相结合提取脑血管图像并分割，包括：S1，数据预处理，得到去噪后脑部区域图像；S2，拟合S1图像背景及血管类的高斯模型，搭建FMM；S3，设置FMM初始值，估计参数及计算低层能量函数；S4，对图像进行分割；S5，结合脊线性质对种子点进行筛选，在3D空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型进行补充。所述方法提高了远端血管图像的分割精度，能实现在较为粗糙的血管图像分割结果下得到精准的血管中心线提取。提取。提取。

【技术实现步骤摘要】
一种基于统计学模型和多尺度滤波的脑血管图像分割方法


[0001]本专利技术涉及一种基于统计学模型和多尺度滤波的脑血管图像分割方法，属于医学计算及图像处理


技术介绍

[0002]血管狭窄、动脉瘤和血管畸形等脑血管疾病是威胁人类健康的重大疾病。随着我国老龄化速度加快，老年脑血管病人数量也在逐年上升，及早发现脑血管形态变化，对脑血管疾病进行有效预防和诊断极为重要。基于脑血管影像的血管数字化分割提取是其中的关键技术。
[0003]脑血管病指的是在大脑中部分血管由于病理原因无法正常运作，使得脑组织缺血而受损的疾病，比较典型的有免疫机制相关的血管炎、肿瘤压迫等导致的血栓、高血压造成的纤维蛋白坏死等，这些疾病大多能够引发病患脑部组织出现缺血或出血，最终导致患者死亡或残疾。据调查，这些疾病的临床患者以中老年人居多，且大多是急性发病。“中风”就是一种脑血管病，临床称之为脑卒中。临床将该疾病划分为两类：出血性卒中与缺血性卒中。调查发现，该类疾病是导致我国成年人残疾甚至死亡的最常见原因。该疾病具有发病率高、复发率高、致残率以及死亡率较高的特点，会对病人家庭造成较重的经济负担，其症状常表现为瘫痪、言语认知出现障碍以及抑郁症等。患病率以及发病率相比于2015年也在持续上升，对于人类的健康造成了巨大的威胁。上文提到的两种脑卒中，治疗方式并不相同。考虑到当前有效治疗方式的缺乏，预防是最为合适的应对措施。由此，对于缺血性卒中等脑血管疾病的相应治疗措施而言，是否能够将完整血管进行三维可视化显得至关重要。
[0004]目前基于深度学习方法进行的血管图像分割已经取得了较好的效果，但其依赖大量的标注数据进行训练，同时需要标准的数据采集规范，限制了模型在临床上的推广应用。而统计分析模型依据像素灰度建立有限混合模型，提取血管过程中无需标注数据，对于小样本的医学数据具有优势，但受其模型限制对于微小血管图像分割效果不理想。本申请致力于解决上述问题。

技术实现思路

[0005]本专利技术的目的是为解决由于依据像素灰度建立有限混合模型的统计分析模型对于高强度区域能有效分割，但对于低强度的小血管图像分割的效果较差这一技术缺陷，提出了一种基于统计学模型和多尺度滤波的脑血管图像分割方法，该方法基于脊线的3D中心线提取，将统计模型和多尺度滤波相结合，通过高层与低层模型相结合实现对脑血管图像提取，具体为：利用统计学模型在高强度像素点上分割提取，结合脊线性质对种子点进行筛选；应用脊线性质，在3D空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，实现了多尺度滤波的脑血管图像分割。
[0006]为了达到上述目的，本专利技术采取如下技术方案。
[0007]所述脑血管图像分割方法，包括如下步骤：
[0008]步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；具体为：
[0009]步骤1.1：通过TOF
‑
MRA及SNAP
‑
MRA采集出大脑部分血管的原始图像；
[0010]其中，TOF
‑
MRA，即飞行时间磁共振血管成像；SNAP
‑
MRA即非对比增强血管与斑块内出血同时成像；血管三维图像即颅骨剥离后的脑部区域图像；
[0011]步骤1.2：使用FSL软件的bet命令对步骤1.1采集的原始图像进行颅骨剥离，将原始数据去除眼球脂肪等不相干的图像，得到血管三维图像；
[0012]其中，步骤1.2中对原始图像进行处理的原因是：TOF
‑
MRA及SNAP
‑
MRA采集时由于设备限制，无法单独采集脑部区域；
[0013]步骤1.3：选用大小为n*n*n，方差为se的高斯核对血管三维图像颅骨剥离后的脑部区域图像进行高斯滤波,滤除图像中的无关噪声，得到去噪后脑部区域图像，n的取值范围是2到10，se的取值范围是0.2
‑
1.0，具体包括如下子步骤：
[0014]步骤1.3.1：构建血管三维图像的黑塞矩阵，其中，血管三维图像中每一个像素点的一阶导数表示灰度的梯度，二阶导数是灰度梯度的梯度，即灰度的曲率；
[0015]步骤1.3.2：求出血管三维图像中每一个像素点的特征向量后，将该像素点的曲率进行分解；判断出当前像素点的曲率主方向及该方向对应的照其绝对值大小顺序排列的特征向量λ1，λ2和λ3；
[0016]步骤1.3.3：建立针对某一像素点(x,y,z)的高斯卷积核G(x,y,z)：
[0017]步骤1.3.4：将高斯卷积核G(x,y,z)与对应像素点(x,y,z)的体素强度 I(x,y,z)进行卷积得到，在尺度σ下的体素强度I
′
(x，y，z，σ)；
[0018]步骤1.3.5重复步骤1.3.3到1.3.4，得到在尺度σ下血管三维图像中所有像素点的体素强度；
[0019]步骤1.3.6：使用Hessian滤波器检测出血管三维图像中的管状结构，再对管状结构进行增强，对其余结构信号进行抑制减弱，以增强血管组织的对比度，得到增强后的血管三维图像，具体为：如果特征值λ2和λ3中都小于0，那么我们引入超参数α用于区分板状结构与管状结构；β用于区分斑点状结构与管状结构；c是用于区分板状结构与斑点状结构，计算对于某一像素点s在尺度σ下增强之后的像素表征V0(s，σ)；
[0020]步骤1.3.7：重复步骤1.3.6，在每一个尺度下σ，进行计算该尺度时各个像素点计算得到的像素表征V0(s，σ)；
[0021]步骤1.3.8：将步骤1.3.7中各个尺度下得到的所有的V0(s，σ)中的最大值作为去噪后血管三维图像该像素点的像素值V(s)；
[0022]其中，去噪后血管三维图像即去噪后脑部区域图像；
[0023]步骤1.3.9：根据步骤1.3.8得到的去噪后的血管三维图像的像素值V(s)的灰度分布情况构建灰度分布直方图I
init
；
[0024]步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型FMM，具体为：
[0025]步骤2.1：使用两个高斯模型G1(x)、G2(x)拟合步骤1.3.8得到的去噪后血管三维图像中的背景类；
[0026]步骤2.2：使用一个高斯模型G3(x)拟合步骤1.3.8得到的去噪后血管三维图像中的血管类；
[0027]步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型G1(x)、G2(x) 以及G3(x)搭建有限混合模型FMM，其FMM公式描述的灰度分布his()是步骤2.1与步骤2.2拟合出来的三个高斯模型的线性组合，his()利用超参数α1，α2，α3调节三个高斯模型G1(x)、G2(x)以及G3(x)间各自的权重；
[0028]步骤3：使用最大期望方法对有本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：包括如下步骤：步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型FMM，具体为：步骤2.1：使用两个高斯模型G1(x)、G2(x)拟合步骤1的去噪后血管三维图像中的背景类；步骤2.2：使用一个高斯模型G3(x)拟合步骤1.3.8得到的去噪后血管三维图像中的血管类；步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型G1(x)、G2(x)以及G3(x)搭建有限混合模型FMM，其FMM公式描述的灰度分布his()是步骤2.1与步骤2.2拟合出来的三个高斯模型的线性组合，his()利用超参数α1，α2，α3调节三个高斯模型G1(x)、G2(x)以及G3(x)间各自的权重；步骤3：使用最大期望方法对有限混合模型FMM设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型FMM设置初始值；步骤4：使用步骤2搭建的有限混合模型FMM与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：步骤4.1：通过步骤3得到的his(y
s
)函数的初始参数θ的估计，提出低层能量函数U1(x
s
)满足U1(x
s
)
∝
his(y
s
)，计算对于去噪后的血管三维图像的单个像素点s的模型G1、G2以及G3的概率表达P(y
s
|G1)、P(y
s
|G2)以及P(y
s
|G3)；步骤4.2：针对去噪后的血管三维图像的单个像素，计算血管类的条件概率函数P(Y|V)和背景类的条件概率函数P(Y|B)；步骤4.3：判断P(Y|V)是否大于P(Y|B)，若是，则认为该像素点属于血管类，否则认为该像素点属于背景类，得到该像素标签；步骤4.4：重复上述步骤4.1
‑
4.3，遍历去噪后血管三维图像的每一个像素，为每一个像素生成其对应的标签；步骤4.5：将步骤4.4中针对每一个像素得到的相应的标签进行汇总，进而得到去噪的血管三维图像分割的初步结果Mask0；步骤5：使用步骤1中去噪后血管三维图像进行基于脊线性质的血管中心线提取，具体为：提取去噪后血管三维图像中的种子点，并且对提取的种子点生成初始方向，之后针对生成的初始方向对每一个种子点沿正向与负向进行延伸以实现中心线生长，具体为：步骤5.1：提取去噪的血管三维图像中的种子点，具体步骤如下：步骤5.1.1：遍历经过步骤1中去噪后血管三维图像中的所有像素点，对于每一个像素点s将该点强度值I
s
与其邻域内的点N(s)强度值进行比较，判断像素点s的灰度是否大于邻域内任意一点灰度，如果是则认为该像素点是一个极值点，否则不认为该像素点是一个极值点；
步骤5.1.2：遍历步骤5.1.1得到的所有极值点，针对每个强度为I0，梯度方向为的极值点s0(x0，y0，z0)进行一阶微分方向修正，找到点s1(x1，y1，z1)，其强度为I1；步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s1(x1，y1，z1)，计算其黑塞矩阵较大的两个特征值对应的特征向量步骤5.1.4：根据步骤5.1.3针对修正过的点s1(x1，y1，z1)及陶成平面Ω，并以s1为圆心r为半径的范围内在平面Ω寻找灰度极值点作为种子点；步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着寻找其邻域内强度最大点s1，判断s1强度是否小于s0，若是则结束此次初始方向的生成，否则将其作为正向传播的下一点s1，称为后沿点；其中，将从种子点s0指向后沿点s1的方向作为正向传播的方向步骤5.2，2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着寻找其邻域内强度最大点s
‑1，判断s
‑1强度小于s0，若是则结束此次初始方向的生成，否则将其作为负向传播的下一点s
‑1，称为前驱点；其中，将从种子点s0指向前驱点s
‑1的方向作为负向传播的方向步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s
‑1从s0的邻域点集内移除，跳至步骤5.2.1；否则完成该种子点的初始方向的生成，跳至步骤5.3；步骤5.3：针对步骤5.2生成的初始方向对每一个种子点沿正向与负向进行延伸实现中心线生长；步骤5.4：对于步骤5.1提取的每一个种子点重复步骤5.2
‑
5.3，得到该血管三维图像的血管中心线图像；步骤5.5：对于步骤5.4中提取到的血管中心线做膨胀处理，得到处理之后的中心线图像；步骤6：引入血管特征场及血管形状场，计算势能函数，并且根据势能函数计算高层能量函数，具体步骤如下：步骤6.1：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一个像素点s，计算该点血管特征场对于高层能量函数的贡献其中，像素点s的邻域点集为N(s)，针对邻域内一点任意一点s
′
∈N(s)使用表示该像素点处梯度，表示其邻域内梯度最大值；步骤6.2：重复步骤6.1，对去噪的血管三维图像数据的每一个像素点求解血管特征场；步骤6.3：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一像素点
s，计算该点血管方向场对于高层能量函数的贡献其中，像素点s的邻域点集为N(s)，使用表示s处梯度，其邻域内一点s
′
∈N(s)，使用表示s
′
处最小特征值对应特征向量，则s
′
处血管方向场对于高层能量函数的贡献为步骤6.4：重复步骤6.3，对去噪的血管三维图像数据的每一个像素点求解血管方向场；步骤6.5：基于步骤6.2和6.4求解的针对所有像素点的血管特征场与血管方向场，对于各个像素点，通过步骤4.5得到的初始标签Mask0计算出该像素点与其邻域内一点s
′
所对应的势能函数Ita(x
s
，x
s
′
)，其为血管特征场与血管方向场的线性组合，用ε1，ε2控制各自的权重；步骤6.6：基于步骤6.5的运算得到的势能函数Ita(x
s
，x
s
′
)，根据结合s
′
处在Mask0中的标签计算对x
s
的能量函数的贡献值E(x
s
，x
s
′
)；步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数U2(x
s
)；步骤7：计算当前分割结果，直到图像分割结果Mask变化小于设立阈值，具体步骤如下：步骤7.1：结合步骤3得到的低层能量函数U1(x
s
)与步骤6.7得到的高层能量函数U2(x
s
)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率P(x
s
)步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次...

【专利技术属性】
技术研发人员：李雪松，姚宇翔，石根，
申请(专利权)人：北京理工大学，
类型：发明
国别省市：