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

权利要求书

1.一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：包括如下步骤：

步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；

步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型FMM，具体为：

步骤2.1：使用两个高斯模型G₁(x)、G₂(x)拟合步骤1的去噪后血管三维图像中的背景类；

步骤2.2：使用一个高斯模型G₃(x)拟合步骤1.3.8得到的去噪后血管三维图像中的血管类；

步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型G₁(x)、G₂(x)以及G₃(x)搭建有限混合模型FMM，其FMM公式描述的灰度分布his()是步骤2.1与步骤2.2拟合出来的三个高斯模型的线性组合，his()利用超参数α₁，α₂，α₃调节三个高斯模型G₁(x)、G₂(x)以及G₃(x)间各自的权重；

步骤3：使用最大期望方法对有限混合模型FMM设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型FMM设置初始值；

步骤4：使用步骤2搭建的有限混合模型FMM与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：

步骤4.1：通过步骤3得到的his(y_s)函数的初始参数θ的估计，提出低层能量函数U₁(x_s)满足U₁(x_s)∝his(y_s)，计算对于去噪后的血管三维图像的单个像素点s的模型G₁、G₂以及G₃的概率表达P(y_s|G₁)、P(y_s|G₂)以及P(y_s|G₃)；

步骤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得到的所有极值点，针对每个强度为I₀，梯度方向为的极值点s₀(x₀，y₀，z₀)进行一阶微分方向修正，找到点s₁(x₁，y₁，z₁)，其强度为I₁；

步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s₁(x₁，y₁，z₁)，计算其黑塞矩阵较大的两个特征值对应的特征向量

步骤5.1.4：根据步骤5.1.3针对修正过的点s₁(x₁，y₁，z₁)及陶成平面Ω，并以s₁为圆心r为半径的范围内在平面Ω寻找灰度极值点作为种子点；

步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：

步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s₀沿着寻找其邻域内强度最大点s₁，判断s₁强度是否小于s₀，若是则结束此次初始方向的生成，否则将其作为正向传播的下一点s₁，称为后沿点；

其中，将从种子点s₀指向后沿点s₁的方向作为正向传播的方向

步骤5.2，2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s₀沿着寻找其邻域内强度最大点s_-1，判断s_-1强度小于s₀，若是则结束此次初始方向的生成，否则将其作为负向传播的下一点s_-1，称为前驱点；

其中，将从种子点s₀指向前驱点s_-1的方向作为负向传播的方向

步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s_-1从s₀的邻域点集内移除，跳至步骤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′)，其为血管特征场与血管方向场的线性组合，用ε₁，ε₂控制各自的权重；

步骤6.6：基于步骤6.5的运算得到的势能函数Ita(x_s，x_s′)，根据结合s′处在Mask0中的标签计算对x_s的能量函数的贡献值E(x_s，x_s′)；

步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数U₂(x_s)；

步骤7：计算当前分割结果，直到图像分割结果Mask变化小于设立阈值，具体步骤如下：

步骤7.1：结合步骤3得到的低层能量函数U₁(x_s)与步骤6.7得到的高层能量函数U₂(x_s)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率P(x_s)

步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次去噪后血管三维图像的分割结果间的相似性系数是否低于阈值c，若是，则输出分割的标签Mask0作为最终的结果，否则Mask0设置为此次去噪的血管三维图像分割的结果并且跳转到对步骤5；

其中，阈值c的取值范围为0.01到0.05。