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

1.本发明涉及一种基于统计学模型和多尺度滤波的脑血管图像分割方法，属于医学计算及图像处理技术领域。


背景技术：

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


技术实现要素：

5.本发明的目的是为解决由于依据像素灰度建立有限混合模型的统计分析模型对于高强度区域能有效分割，但对于低强度的小血管图像分割的效果较差这一技术缺陷，提出了一种基于统计学模型和多尺度滤波的脑血管图像分割方法，该方法基于脊线的3d中心线提取，将统计模型和多尺度滤波相结合，通过高层与低层模型相结合实现对脑血管图像提取，具体为：利用统计学模型在高强度像素点上分割提取，结合脊线性质对种子点进行筛选；应用脊线性质，在3d空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，实现了多尺度滤波的脑血管图像分割。
6.为了达到上述目的，本发明采取如下技术方案。
7.所述脑血管图像分割方法，包括如下步骤：
8.步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；具体为：
9.步骤1.1：通过tof-mra及snap-mra采集出大脑部分血管的原始图像；
10.其中，tof-mra，即飞行时间磁共振血管成像；snap-mra即非对比增强血管与斑块内出血同时成像；血管三维图像即颅骨剥离后的脑部区域图像；
11.步骤1.2：使用fsl软件的bet命令对步骤1.1采集的原始图像进行颅骨剥离，将原始数据去除眼球脂肪等不相干的图像，得到血管三维图像；
12.其中，步骤1.2中对原始图像进行处理的原因是：tof-mra及snap-mra采集时由于设备限制，无法单独采集脑部区域；
13.步骤1.3：选用大小为n*n*n，方差为se的高斯核对血管三维图像颅骨剥离后的脑部区域图像进行高斯滤波,滤除图像中的无关噪声，得到去噪后脑部区域图像，n的取值范围是2到10，se的取值范围是0.2-1.0，具体包括如下子步骤：
14.步骤1.3.1：构建血管三维图像的黑塞矩阵，其中，血管三维图像中每一个像素点的一阶导数表示灰度的梯度，二阶导数是灰度梯度的梯度，即灰度的曲率；
15.步骤1.3.2：求出血管三维图像中每一个像素点的特征向量后，将该像素点的曲率进行分解；判断出当前像素点的曲率主方向及该方向对应的照其绝对值大小顺序排列的特征向量λ1，λ2和λ3；
16.步骤1.3.3：建立针对某一像素点(x,y,z)的高斯卷积核g(x,y,z)：
17.步骤1.3.4：将高斯卷积核g(x,y,z)与对应像素点(x,y,z)的体素强度 i(x,y,z)进行卷积得到，在尺度σ下的体素强度i
′
(x，y，z，σ)；
18.步骤1.3.5重复步骤1.3.3到1.3.4，得到在尺度σ下血管三维图像中所有像素点的体素强度；
19.步骤1.3.6：使用hessian滤波器检测出血管三维图像中的管状结构，再对管状结构进行增强，对其余结构信号进行抑制减弱，以增强血管组织的对比度，得到增强后的血管三维图像，具体为：如果特征值λ2和λ3中都小于0，那么我们引入超参数α用于区分板状结构与管状结构；β用于区分斑点状结构与管状结构；c是用于区分板状结构与斑点状结构，计算对于某一像素点s在尺度σ下增强之后的像素表征v0(s，σ)；
20.步骤1.3.7：重复步骤1.3.6，在每一个尺度下σ，进行计算该尺度时各个像素点计算得到的像素表征v0(s，σ)；
21.步骤1.3.8：将步骤1.3.7中各个尺度下得到的所有的v0(s，σ)中的最大值作为去噪后血管三维图像该像素点的像素值v(s)；
22.其中，去噪后血管三维图像即去噪后脑部区域图像；
23.步骤1.3.9：根据步骤1.3.8得到的去噪后的血管三维图像的像素值v(s)的灰度分布情况构建灰度分布直方图i
init
；
24.步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型fmm，具体为：
25.步骤2.1：使用两个高斯模型g1(x)、g2(x)拟合步骤1.3.8得到的去噪后血管三维图像中的背景类；
26.步骤2.2：使用一个高斯模型g3(x)拟合步骤1.3.8得到的去噪后血管三维图像中的血管类；
27.步骤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)间各自的权重；
28.步骤3：使用最大期望方法对有限混合模型fmm设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型fmm设置初始值，具体为：
29.步骤3.1：为有限混合模型设置初始值，得到三个高斯模型的参数集θ，具体包括如下子步骤：
30.步骤3.1.1：初始化步骤2.2拟合得到的血管类的高斯分布g3(x)，将其在步骤2.3中得到的his(x)表达式中的占比α3设定为区间0.01～0.05中的常数值；
31.步骤3.1.2：使用经过步骤1.3.8中去噪的血管三维图像中最亮的前m％的像素点计算得出经过步骤2.2拟合出来的高斯分布函数g3(x)的初始均值μ3及方差σ
32
；
32.其中，m为取值1-5之间的常数；
33.步骤3.1.3：将经过步骤2.1拟合出来的g2(x)的初始均值μ2设定为步骤 1.3.9中得到的灰度直方图i
init
的最高点的灰度；
34.步骤3.1.4：将经过步骤2.1拟合出来的g1(x)的初始均值μ1设立为g2(x)的初始均值μ2的t倍；
35.其中，t为常数，取值范围为0.1-0.5；
36.步骤3.1.5：计算g1(x)贡献在步骤1.3.9中的得到的灰度直方图中i
init
占比作为α1的估计值；
37.步骤3.1.6：根据步骤3.1.5计算得到的g1(x)的体素占比α1，计算中去掉 g1(x)贡献后的剩余分布直方图i
left
；
38.步骤3.1.7：取步骤3.1.6得到的剩余分布直方图i
left
在区间[(μ2-μ1)/ 2，(μ2+μ1)/2]上的灰度点进行最大似然估计计算g2(x)函数的方差σ
22
；
[0039]
步骤3.1.8：根据步骤3.1.1的求解的α3和步骤3.1.5的求解的α1，g2(x) 的占比α2使用1-α
1-α3得到；
[0040]
均值参数的集合θ
μ
、方差参数的的集合θ
σ
及体素占比参数的集合θ
α
，统称为参数集θ；
[0041]
步骤3.2：对有限混合模型执行最大期望方法，具体步骤如下：
[0042]
步骤3.2.1：设定循环参数t＝1；
[0043]
步骤3.2.2：根据三个高斯模型的参数集θ，计算隐变量γ其联合概率的条件概率p；
[0044]
步骤3.2.3：根据步骤3.2.2计算出来的隐变量γ其联合概率的条件概率p，计算期望e；
[0045]
步骤3.2.4：根据步骤3.2.3计算得到的期望e，计算使得其最大的参数取值，作为t+1次循环的参数集，将体素占比的和，即∑θ
α
设为1，得到新拉格朗日函数；
[0046]
步骤3.2.5：使用步骤3.2.4得到的新拉格朗日函数对λ、αz、μz、σz求偏导，令其偏导数为0，得到新的参数集；
[0047]
步骤3.2.6：判断t是否大于等于2，若是，则判断第t次及第t-1两次循环的所有参数的变化是否均小于q％，若是，将步骤3.2.5得到的新的参数集作为 his(ys)函数的初始参数的估计，并跳至步骤4；否则判断t是否等于n，若是，将步骤3.2.5得到的新的参数集作为his(ys)函数的初始参数的估计，并跳至步骤4，否则将t加1，跳至步骤3.2.2；
[0048]
其中，q的取值范围在0.5-1.5之间，n的取值范围为25到75；
[0049]
至此，完成了对his(ys)函数的初始参数θ的估计；
[0050]
步骤4：使用步骤2搭建的有限混合模型fmm与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：
[0051]
步骤4.1：通过步骤3得到的his(ys)函数的初始参数θ的估计，提出低层能量函数u1(xs)满足u1(xs)
∝
his(ys)，计算对于去噪后的血管三维图像的单个像素点s的模型g1、g2以及g3的概率表达p(ys|g1)、p(ys|g2)以及p(ys|g3)；
[0052]
步骤4.2：针对去噪后的血管三维图像的单个像素，计算血管类的条件概率函数p(y|v)和背景类的条件概率函数p(y|b)；
[0053]
步骤4.3：判断p(y|v)是否大于p(y|b)，若是，则认为该像素点属于血管类，否则认为该像素点属于背景类，得到该像素标签；
[0054]
步骤4.4：重复上述步骤4.1-4.3，遍历去噪后血管三维图像的每一个像素，为每一个像素生成其对应的标签；
[0055]
步骤4.5：将步骤4.4中针对每一个像素得到的相应的标签进行汇总，进而得到去噪的血管三维图像分割的初步结果mask0；
[0056]
步骤5：使用步骤1中去噪后血管三维图像进行基于脊线性质的血管中心线提取，具体为：提取去噪后血管三维图像中的种子点，并且对提取的种子点生成初始方向，之后针对生成的初始方向对每一个种子点沿正向与负向进行延伸以实现中心线生长，具体为：
[0057]
步骤5.1：提取去噪的血管三维图像中的种子点,具体步骤如下：
[0058]
步骤5.1.1：遍历经过步骤1中去噪后血管三维图像中的所有像素点，对于每一个像素点s将该点强度值is与其邻域内的点n(s)强度值进行比较，判断像素点s的灰度是否大于邻域内任意一点灰度，如果是则认为该像素点是一个极值点，否则不认为该像素点是一个极值点；
[0059]
步骤5.1.2：遍历步骤5.1.1得到的所有极值点，针对每个强度为i0，梯度方向为的极值点s0(x0，y0，z0)进行一阶微分方向修正，找到点 s1(x1，y1，z1)，其强度为i1；
[0060]
步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s1(x1，y1，z1)，计算其黑塞矩阵较大的两个特征值对应的特征向量
[0061]
步骤5.1.4：根据步骤5.1.3针对修正过的点s1(x1，y1，z1)及构成平面ω，并以s1为圆心r为半径的范围内在平面ω寻找灰度极值点作为种子点；
[0062]
步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：
[0063]
步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点 s0沿着寻找其邻域内强度最大点s1，判断s1强度是否小于s0，若是则结束此次初始方向的生成，否则将其作为正向传播的下一点s1，称为后沿点；
[0064]
其中，将从种子点s0指向后沿点s1的方向作为正向传播的方向
[0065]
步骤5.2,2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着寻找其邻域内强度最大点s-1
，判断s-1
强度小于s0，若是则结束此次初始方向的生成，否则将其作为负向传播的下一点s-1
，称为前驱点；
[0066]
其中，将从种子点s0指向前驱点s-1
的方向作为负向传播的方向
[0067]
步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s-1
从s0的邻域点集内移除，跳至步骤5.2.1；否则完成该种子点的初始方向的生成，跳至步骤5.3；
[0068]
步骤5.3：针对步骤5.2生成的初始方向对每一个种子点沿正向与负向进行延伸实现中心线生长；
[0069]
步骤5.4：对于步骤5.1提取的每一个种子点重复步骤5.2—5.3，得到该血管三维图像的血管中心线图像；
[0070]
步骤5.5：对于步骤5.4中提取到的血管中心线做膨胀处理，得到处理之后的中心线图像；
[0071]
步骤6：引入血管特征场及血管形状场，计算势能函数，并且根据势能函数计算高层能量函数，具体步骤如下：
[0072]
步骤6.1：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一个像素点s，计算该点血管特征场对于高层能量函数的贡献
[0073]
其中，像素点s的邻域点集为n(s)，针对邻域内一点任意一点s
′
∈n(s)使用表示该像素点处梯度，表示其邻域内梯度最大值；
[0074]
步骤6.2：重复步骤7.1，对去噪的血管三维图像数据的每一个像素点求解血管特征场；
[0075]
步骤6.3：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一像素点s，计算该点血管方向场对于高层能量函数的贡献
[0076]
其中，像素点s的邻域点集为n(s)，使用表示s处梯度，其邻域内一点s
′
∈ n(s)，使用表示s
′
处最小特征值对应特征向量，则s
′
处血管方向场对于高层能量函数的贡献为
[0077]
步骤6.4：重复步骤6.3，对去噪的血管三维图像数据的每一个像素点求解血管方向场；
[0078]
步骤6.5：基于步骤6.2和6.4求解的针对所有像素点的血管特征场与血管方向场，对于各个像素点，通过步骤4.5得到的初始标签mask0计算出该像素点与其邻域内一点s
′
所对应的势能函数ita(xs,xs′
)，其为血管特征场与血管方向场的线性组合，用ε1,ε2控制各自的权重；
[0079]
步骤6.6：基于步骤6.5的运算得到的势能函数ita(xs,xs′
)，根据结合s
′
处在 mask0中的标签计算对xs的能量函数的贡献值e(xs,x
s'
)；
[0080]
步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数u2(xs)；
[0081]
步骤7：计算当前分割结果，直到图像分割结果mask变化小于设立阈值，具体步骤如下:
[0082]
步骤7.1：结合步骤3得到的低层能量函数u1(xs)与步骤6.7得到的高层能量函数u2(xs)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率 p(xs)
[0083]
步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次去噪后血管三维图像的分割结果间的相似性系数是否低于阈值c，若是，则输出分割的标签 mask0作为最终的结果，否则mask0设置为此次去噪的血管三维图像分割的结果并且跳转到对步骤5；
[0084]
其中，阈值c的取值范围为0.01到0.05。
[0085]
有益效果
[0086]
本发明一种基于统计学模型和多尺度滤波的tof-mra脑血管图像分割方法，与已有技术相比较，具有如下有益效果：
[0087]
1.所述血管图像分割方法针对统计学模型中应用脊线性质，在3d空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，提高低层模型对远端血管图像分割精度，解决了低层模型对于远端血管图像分割效果较差的问题；
[0088]
2.所述方法基于脊线的3d中心线提取，即利用统计学模型在高强度像素点上分割提取的优势，结合脊线性质对种子点进行筛选，可以实现在较为粗糙的血管图像分割结果下得到精准的血管中心线提取。
附图说明
[0089]
图1为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法的中心线膨胀示意图；
[0090]
图2为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法的流程图；
[0091]
图3为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法中地面真实标注图像的示意图；
[0092]
图4为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法一张切片下原图像(a)与增强图像(b)的对比示意图；
[0093]
图5中应用本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法给出了颅骨剥离后的血管图像灰度分布直方图i
init
；
[0094]
图6中应用本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法低层模型结果对比；
[0095]
图7中应用本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法种子点延伸示意图对比；
[0096]
图8为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法中心线分割的两例结果对比图；
[0097]
图9为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法针对brave数据集的运行结果。
具体实施方式
[0098]
下面结合附图和实例，对本发明提供的一种基于统计学模型和多尺度滤波的脑血管图像分割方法作详细地说明。
[0099]
实施例1
[0100]
脑血管疾病是导致人类死亡的主要疾病之一，对其早期预防及治疗极为重要。通过影像技术进行的无创伤脑血管病变的检测，是一种有效的诊断和监测脑血管疾病的方法，也是最能被人所接受而不会产生副作用的技术。此项研究为脑血管疾病的诊断与预防提供技术支持，相关技术可以推广到多模式图像配准和三维可视化等其它领域。由于我国每年参加体检和预防性保健人数规模庞大，此研究成果的应用可产生广泛的社会效益和经济效益。
[0101]
本实施例中以brave数据集作为我们的目标对象；brave，即brain andvascular health in the elderly，此数据集研究的目的是使用多模式核磁功能成像来研究血管健康和大脑老化之间的相互关系。这项研究连续招募了200名老年男性，他们患有原发性高血压，没有临床上明显的神经、心脑血管疾病。所有与会者都提供了书面知情同意。
[0102]
数据集在进行采集时，由于设备原因使用了分段扫描并拼接，使得整张图象在成像时呈现了明显的灰度阶梯，这部分噪声使得最终成像效果较差，进而影响到了分割结果，使得整体参数较低。由于医疗数据获取地面真实困难，对于一例 tof-mra的血管进行标注就需要一位有经验的外科医生3-4天的工作量，本课题所用数据集对于血管分割部分仅标注了关心的血管结果，对于其余血管并未进行标注。具体标注结果如附图3所示，在附图3中左侧图为tof-mra切片，右侧图为对应切片的地面真实与切片叠加展示，使用红色表示地面真实的标注，本实施例在brave数据集中随机选取了20例数据进行试验；结合实际情况，对于血管分割结果本课题使用了真实阳性值作为核心指标。
[0103]
地面真实，也即ground truth，即为真值、真实的有效值或者是标准的答案。
[0104]
本实施例中的一种基于统计学模型和多尺度滤波的tof-mra脑血管图像分割方法，其操作流程如附图2所示，具体实现步骤为：
[0105][0106]
步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；具体为：
[0107]
步骤1.1：通过tof-mra及snap-mra采集出大脑部分血管的原始图像；
[0108]
其中，tof-mra，即飞行时间磁共振血管成像；snap-mra即非对比增强血管与斑块内出血同时成像；血管三维图像即颅骨剥离后的脑部区域图像；
[0109]
步骤1.2：使用fsl软件的bet命令对步骤1.1采集的原始图像进行颅骨剥离，将原始数据去除眼球脂肪等不相干的图像，得到血管三维图像；
[0110]
其中，步骤1.2中对原始图像进行处理的原因是：tof-mra及snap-mra采集时由于设备限制，无法单独采集脑部区域；
[0111]
步骤1.3：选用大小为3
×3×
3，方差为0.4的高斯核对血管三维图像颅骨剥离后的脑部区域图像进行高斯滤波,滤除图像中的无关噪声，得到去噪后脑部区域图像，具体包括如下子步骤：
[0112]
步骤1.3.1：构建血管三维图像的黑塞矩阵，其中，血管三维图像中每一个像素点的一阶导数表示灰度的梯度，二阶导数是灰度梯度的梯度，即灰度的曲率，其黑塞矩阵的表达形式为公式(1)：
[0113][0114]
步骤1.3.2：求出血管三维图像中每一个像素点的特征向量后，将该像素点的曲率进行分解；判断出当前像素点的曲率主方向及该方向对应的照其绝对值大小顺序排列的特征向量λ1，λ2和λ3；
[0115]
步骤1.3.3：建立针对某一像素点(x,y,z)的高斯卷积核g(x,y,z)，如公式 (2)所示：
[0116][0117]
步骤1.3.4将高斯卷积核g(x,y,z)与对应像素点(x,y,z)的体素强度 i(x,y,z)进行卷积得到，在尺度σ下的体素强度i
′
(x，y，z，σ)；
[0118]
步骤1.3.5重复步骤1.3.3到1.3.4，得到在尺度σ下血管三维图像中所有像素点的体素强度，计算方法如公式(3)所示：
[0119]i′
(x,y,z,σ)＝g(x,y,z,σ)*i(x,y,z)
ꢀꢀꢀ
(3)
[0120]
步骤1.3.6：使用hessian滤波器检测出血管三维图像中的管状结构，再对管状结构进行增强，对其余结构信号进行抑制减弱，以增强血管组织的对比度，得到增强后的血管三维图像，具体为：如果特征值λ2和λ3中都小于0，那么我们引入超参数α用于区分板状结构与管状结构；β用于区分斑点状结构与管状结构；c是用于区分板状结构与斑点状结构，计算对于某一像素点s在尺度σ下增强之后的像素表征v0(s，σ)，如公式(4)所示：
[0121]
[0122]
步骤1.3.7：重复步骤1.3.6，在每一个尺度下σ，进行计算该尺度时各个像素点计算得到的结果；
[0123]
步骤1.3.8：将步骤1.3.7中各个尺度下得到的所有的v0(s，σ)中的最大值作为去噪后血管三维图像该像素点的像素值v(s)，如公式(5)所示：
[0124][0125]
其中，去噪后血管三维图像即去噪后脑部区域图像；
[0126][0127][0128]
附图4给出了图像中一张切片下原图像(a)与增强图像(b)的对比。图像经过增强后，可以鲜明的看到其中的血管部分被增强，但其中也有少量的非血管组织也被错误的增强；同时由于血管尺度不同，其增强的血管会有部分变得比原图像更加细，而由于这样的原因多尺度滤波器更多的是作为血管图像分割时图像的预处理手段而不是直接作为分割结果；
[0129]
步骤1.3.9：根据步骤1.3.8得到的处理之后的血管三维图像的像素图v(s) 的灰度分布情况构建灰度分布直方图i
init
；
[0130][0131]
附图5中给出了颅骨剥离后的血管图像灰度分布直方图i
init
：纵坐标表示频率，横坐标表示图像像素的灰度；在灰度分布直方图中，最低强度级别主要对应于脑脊液、骨骼和背景空气；中等强度级别对应于脑组织，包括灰质和白质以及眼睛的部分；高强度级别由邻近皮肤的脂肪和动脉组成；基于tof-mra的成像原理，血液流速越快的部分强度越高，即大血管理论上是图像中最亮的部分，与周围组织对比度较高，也更加容易通过灰度的差异来进行分离；在血管末梢部位血管变细，血液流速降低，对应部位在图像中的强度、与周围组织的对比度也会下降，较难通过灰度差异分离；在低层图像处理时，使用的是颅骨剥离后的图像，这是由于骨骼部分对应的强度范围内几乎不会有血管体素的存在，然而低强度区域的灰度分布在拟合时会引入新的误差，进而影响中强度区域及高强度区域的拟合效果；同时去掉颅骨区域后，减少了灰度拟合时的模型个数，更容易拟合；
[0132]
步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型fmm，具体为：
[0133]
步骤2.1：使用两个高斯模型g1(x)、g2(x)拟合步骤1.3.8得到的去噪后血管三维图像中的背景类的背景类，其中α1、μ1、σ
12
代表g1(x)中的权重、均值与方差，α2、μ2、σ
22
代表g2(x)中的权重、均值与方差。对于一个高斯模型gn(x)，其某像素点x(i)处的公式表达为公式(6)：
[0134][0135]
步骤2.2：单独使用一个高斯模型g3(x)拟合其血管类，其参数表示为α3、μ3、σ
32
，其高斯函数的表达形式如公式(6)所示。
[0136]
步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型g1(x)、g2(x) 以及g3(x)搭建有限混合模型fmm以拟合灰度分布his()，公式如(7)所示：
[0137]
his(x)＝α1×
g1(x)+α2×
g2(x)+α3×
g3(x)
ꢀꢀ
(7)
[0138]
其中超参数α1，α2，α3用于调节步骤2.1与步骤2.2拟合出来的三个高斯模型之间的权重；
[0139]
步骤3：使用最大期望方法对有限混合模型fmm设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型fmm设置初始值，具体为：
[0140]
步骤3.1：为有限混合模型设置初始值，得到三个高斯模型的参数集θ，具体包括如下子步骤：
[0141]
步骤3.1.1：首先初始化步骤2.2拟合得到的血管类的高斯分布g3(x)，将其在步骤2.3中得到的his(x)表达式中的占比α3设定为区间0.03；
[0142]
其中，步骤3.1.1中将超参数α3设定为0.03的原因是：人的大脑中血管体素所占比例通常为0.01～0.05，将其设立为0.03可以较好的进行拟合；
[0143]
步骤3.1.2：将经过步骤2.2拟合出来的高斯分布函数g3(x)的初始均值μ3及方差σ
32
使用经过步骤1.3.8的中经去噪的血管三维图像中的最亮的前3％点进行计算得出；
[0144]
步骤3.1.3：将经过步骤2.1拟合出来的g2(x)的初始均值μ2设定为步骤 1.3.9中得到的灰度直方图i
init
的最高点的灰度；
[0145]
步骤3.1.4：将经过步骤2.1拟合出来的g1(x)的初始均值μ1设立为g2(x)的初始均值μ2的0.3倍；
[0146]
步骤3.1.5：计算g1(x)贡献在步骤1.3.9中的得到的灰度直方图中i
init
占比作为α1的估计值，引入超参数k制衡g1(x)的峰值高度，计算方法参照公式(8)：
[0147][0148]
步骤3.1.6：根据步骤3.1.5计算得到的结果，计算中去掉g1(x)贡献后的剩余分布直方图i
left
，计算方法参照公式(9)：
[0149][0150]
步骤3.1.7：取步骤3.1.6得到的剩余分布直方图i
left
在区间[(μ2-μ1)/ 2，(μ2+μ1)/2]上的灰度点进行最大似然估计计算方差σ
22
；
[0151]
步骤3.1.8：根据步骤3.1.1的求解的α3和步骤3.1.5的求解的α1，g2(x) 的占比α2使用1-α
1-α3得到；
[0152]
均值参数的集合θ
μ
、方差参数的的集合θ
σ
及体素占比参数的集合θ
α
，统称为参数集θ；
[0153]
步骤3.2：对有限混合模型执行最大期望方法，具体步骤如下：
[0154]
步骤3.2.1：设定循环参数t＝1；
[0155]
步骤3.2.2：根据三个高斯模型的参数集θ，计算隐变量γ其联合概率的条件概率p，计算方法参照公式(10)；
[0156][0157]
步骤3.2.3：根据步骤3.2.2计算出来的隐变量γ其联合概率的条件概率p，计算期望e，如公式(11)所示：
[0158][0159]
步骤3.2.4：根据步骤3.2.3计算得到的期望e，计算使得其最大的参数取值，作为t+1次循环的参数集，将体素占比的和，即∑θ
α
设为1，得到新拉格朗日函数l，如公式(12)所示：
[0160][0161]
步骤3.2.4：使用步骤3.2.3得到的l对λ、αz、μz、σz求偏导，令其偏导数为0，可得新的参数集θ
′
，如公式(13)所示；
[0162][0163]
步骤3.2.5：使用步骤3.2.4得到的新拉格朗日函数对λ、αz、μz、σz求偏导，令其偏导数为0，得到新的参数集；
[0164]
步骤3.2.6：判断t是否大于等于2，若是，则判断第t次及第t-1两次循环的所有参数的变化是否均小于1％，若是，将步骤3.2.5得到的新的参数集作为 his(ys)函数的初始参数的估计，并跳至步骤4；否则判断t是否等于50，若是，将步骤3.2.5得到的新的参数集作为his(ys)函数的初始参数的估计，并跳至步骤4，否则将t加1，跳至步骤3.2.2；
[0165]
至此，完成了对his(ys)函数的初始参数θ的估计；
[0166]
步骤4：使用步骤2搭建的有限混合模型fmm与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：
[0167]
步骤4.1：通过步骤3得到的his(ys)函数的初始参数θ的估计，提出低层能量函数u1(xs)满足u1(xs)
∝
his(ys)，计算对于去噪后的血管三维图像的单个像素点s的模型g1、g2以及g3的概率表达p(ys|g1)、p(ys|g2)以及p(ys|g3)。对于单个像素点s，计算其属于第i个模型gi的概率表达p(ys|gi)，如公式(14)所示：
[0168][0169]
步骤4.2：针对去噪后的血管三维图像的单个像素，计算血管类的条件概率函数p(y|v)＝p(ys|g3)和背景类的条件概率函数p(y|b)＝p(ys|g1)+p(ys|g2)；
[0170]
步骤4.3：判断p(y|v)是否大于p(y|b)，若是，则认为该像素点属于血管类，否则认
为该像素点属于背景类，得到该像素标签；
[0171]
步骤4.4：重复上述步骤4.1-4.3，遍历去噪后血管三维图像的每一个像素，为每一个像素生成其对应的标签；
[0172]
步骤4.5：将步骤4.4中针对每一个像素得到的相应的标签进行汇总，进而得到去噪的血管三维图像分割的初步结果mask0；
[0173][0174]
附图6给出了两例病人的通过底层模型的血管分割的对比图，每一行代表一例图像，左侧图是三高斯模型的生成结果、中间图为分段三高斯模型结果，右侧图为标准结果；图像结果均使用最大密度投影图呈现；
[0175]
步骤5：使用步骤1中去噪后血管三维图像进行基于脊线性质的血管中心线提取，具体为：提取去噪后血管三维图像中的种子点，并且对提取的种子点生成初始方向，之后针对生成的初始方向对每一个种子点沿正向与负向进行延伸以实现中心线生长，具体为：
[0176]
步骤5.1：提取去噪的血管三维图像中的种子点,具体步骤如下：
[0177]
步骤5.1.1：遍历经过步骤1中去噪后血管三维图像中的所有像素点，对于一个像素点s，将该点强度值is与其邻域内的点n(s)强度值进行比较，判断像素点 s的灰度是否大于邻域内任意一点灰度，如果是则认为该像素点是一个极值点，否则不认为该像素点是一个极值点；
[0178]
步骤5.1.2：遍历步骤5.1.1得到的所有极值点，针对每个强度为i0，梯度方向为的极值点s0(x0，y0，z0)进行一阶微分方向修正，找到点 s1(x1，y1，z1)，其强度为i1，如公式(15)所示；
[0179][0180]
步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s1(x1，y1，z1)，计算其黑塞矩阵较大的两个特征值对应的特征向量
[0181]
步骤5.1.4：根据步骤5.1.3针对修正过的点s1(x1，y1，z1)及构成平面ω，并以s1为圆心r为半径的范围内在平面ω寻找灰度极值点作为种子点，计算方法参考公式(16)；
[0182][0183]
步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：
[0184]
步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点 s0沿着寻找其邻域内强度最大点s1，若s1强度小于s0，则结束此次初始方向的生成，否则将其作为正向传播的下一点s1，称为后沿点；
[0185]
其中，将从种子点s0指向后沿点s1的方向作为正向传播的方向
[0186]
步骤5.2.2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着
寻找其邻域内强度最大点s-1
，若s-1
强度小于s0，则结束此次初始方向的生成，否则将其作为负向传播的下一点s-1
，称为前驱点；
[0187]
其中，将从种子点s0指向前驱点s-1
的方向作为负向传播的方向
[0188]
步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s-1
从s0的邻域点集内移除，跳至步骤5.2.1；否则完成该种子点的初始方向的生成，跳至步骤5.3；
[0189][0190]
种子点初始方向生成过程如附图7所示，其中p0点表示当前点，对当前点有正向传播及负向传播，直到延伸至强度达到生长阈值或是强度开始下降，则停止生长。
[0191]
步骤5.3：针对步骤5.2生成的初始方向对每一个种子点沿正向与负向进行延伸实现中心线生长；
[0192]
其中，步骤5.2-5.3应用了脊线性质，在3d空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，提高低层模型对远端血管图像分割精度，解决了低层模型对于远端血管图像分割效果较差的问题，以达到我们预期的有益效果；
[0193]
步骤5.4：对于步骤5.1提取的每一个种子点重复步骤5.2—5.3，得到该血管三维图像的血管中心线图像；
[0194][0195]
附图8给出了两例结果对比图，其中第一列为本专利所述中心线提取方法，第三列为中心线的标准分割结果，每一行代表了一例不同的图像。
[0196]
步骤5.5：对于步骤5.4中提取到的血管中心线做膨胀处理，得到处理之后的中心线图像，对于步骤5.4中提取到的血管中心线与血管图像分割的结果相交点，如附图1中点p，若其在实际图像中对应点p0处半径为r，则对中心线pq 做膨胀处理，使得血管半径沿pq逐渐变小，直到汇聚于pq延长线。若pq长度为d，对于pq中与p距离为d的一点，其则其膨胀核尺寸r满足因中心线点不会出现在血管末梢，因此对于r加上限制r＝max(r,2)；
[0197][0198]
步骤6：引入血管特征场及血管形状场，计算势能函数，并且根据势能函数计算高层能量函数，具体步骤如下：
[0199]
步骤6.1：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一个像素点s，计算该点血管特征场对于高层能量函数的贡献
[0200]
其中，像素点s的邻域点集为n(s)，针对邻域内一点任意一点s
′
∈n(s)使用表示该像素点处梯度，表示其邻域内梯度最大值；
[0201]
步骤6.2：重复步骤7.1，对去噪的血管三维图像数据的每一个像素点求解血管特征场；
[0202]
步骤6.3：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一像
素点s，计算该点血管方向场对于高层能量函数的贡献
[0203]
其中，像素点s的邻域点集为n(s)，使用表示s处梯度，其邻域内一点s
′
∈ n(s)，使用表示s
′
处最小特征值对应特征向量，则s
′
处血管方向场对于高层能量函数的贡献为
[0204]
步骤6.4：重复步骤6.3，对去噪的血管三维图像数据的每一个像素点求解血管方向场；
[0205]
步骤6.5：基于步骤6.2和6.4求解的针对所有像素点的血管特征场与血管方向场，对于各个像素点，通过步骤4.5得到的初始标签mask0计算出该像素点与其邻域内一点s
′
所对应的势能函数ita(xs,xs′
)，其为血管特征场与血管方向场的线性组合，用ε1,ε2控制各自的权重，如公式(17)所示：
[0206][0207]
步骤6.6：基于步骤6.5的运算得到的势能函数ita(xs,xs′
)，根据结合s
′
处在 mask0中的标签计算对xs的能量函数的贡献值e(xs,x
s'
)，如公式(18)所示：
[0208][0209]
步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数u2(xs)，通过设置常数βc其随着团结构的不同及图像质量不同而变化。对图像内所有点进行计算，分别记录其对应lv以及lb的能量函数，并将其用于最终血管图像分割结果的计算，如公式(19)所示；
[0210][0211]
步骤7：计算当前分割结果，直到图像分割结果mask变化小于设立阈值，具体步骤如下:
[0212]
步骤7.1：结合步骤3得到的低层能量函数u1(xs)与步骤6.7得到的高层能量函数u2(xs)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率 p(xs)，如公式(20)所示:
[0213][0214]
步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次去噪后血管三维图像的分割结果间的相似性系数是否低于阈值0.03，若是，则输出分割的标签mask0作为最终的结果，否则mask0设置为此次去噪的血管三维图像分割的结果并且跳转到对步骤5；
[0215]
针对brave数据集，上述方法的运行结果如附图9所示：
[0216][0217]
图9中的brave数据集实验结果展示了对比第一列为本课题所提方法结果，第二列血管提取的标准结果。分割结果图像均使用最大密度投影图给出。针对 brave数据集对比结果，可以发现针对图像质量较低、信噪比较差图像，本专利所提方法对血管图像分割准确性有实际提升，尤其对于统计学模型表现不佳的细小血管部分，由于中心线方法的引入，使得该部分结果相较于传统统计学模型有较大的改善，细小血管的分割结果变得更粗，证明了其有效性。
[0218]
上述描述对本发明的特征和方法进行了具体的说明，但应了解，在所述权利要求中定义的本发明并不局限于所述的具体特征或方法。本领域人员可在权利要求的范围内做出修改，并不影响本发明的实质内容。