具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

参照图1，本发明提供了一种脑小血管结构特征的提取方法，包括：

S1：获取脑血管图像，并对获取的图像进行预处理，得到脑小血管区域图像：

本实施例，通过核磁共振成像MRI获取脑血管图像；本实施例所用MRI为7T MRI；

采用linux系统虚拟机中的FSL5.0软件对获取的图像进行颅骨剥离(BET brainextraction)；

对颅骨剥离后的图像进行非均匀校正，其中非均匀校正采用Python中的SimpleITK包进行N4偏置场校正；

对非均匀校正后的图像进行0-255标准化处理，得到处理后图像；

对处理后图像进行区域选定，选定的区域为双侧大脑中动脉MCA区域，该区域包括大脑中动脉MCA起始处到分叉处区域和从MCA发出的所有脑小血管起始到终末区域；

采用灰度阈值法分割选定的区域(即双侧MCA区域)；其中设置血管下限、背景上限两个阈值，基于血管下限和背景上限两个阈值将选定区域分割为背景区域、大血管区域、脑小血管区域；其中灰度值小于背景上限则为背景区域，灰度值大于血管下限则为大血管区域，灰度值在背景上限和血管下限之间则为脑小血管区域；

图像中还存在无法被分割的散在血管点，根据聚类分析cluster算法和势能算法来确定散在血管点，并将确定的散在血管点纳入到分割出来的脑小血管区域中；其中确定散在血管点的方式为：根据连通性相关的cluster算法，得到多个分离的数据簇cluster，以最大的cluster为基准，在剩余的cluster中去除掉体素(voxel)小于100的和Hessian矩阵三个特征值中不符合管状结构的；其中通过Hessian矩阵的三个特征值λ₁、λ₂、λ₃判断该点是否满足管状结构，当Hessian矩阵的特征值λ₁/λ₂<20，返回λ₃的绝对值abs(λ₃)>0.05时，该点不符合管状结构；且设定Hessian矩阵的大小为15×15×15，且三个特征值的绝对值按大小排列后λ₁＞＞λ₂，λ₃≈0；随后采用势能算法，遍历最大的cluster周围的点，将散在血管信号纳入到分割出来的血管区域当中，如果该点满足：在周围3×3×3的区域内亮度大于周围26个点平均值以及至少16个点的一定倍数(本实施例中为确保局部灰度值突出，设置16个点的1.5-2倍，其余图像此处倍数应视其图像质量而定)。且当该点被设定为血管点后，Hessian矩阵仍满足，则认定是血管点；迭代上述过程，迭代次数视图像血管连接度情况而定；迭代步骤，在健康人中可以跳过；

基于上述步骤，可获取到完整的脑小血管区域图像。

基于完整的脑小血管区域生成血管掩膜，生成血管掩膜后，利用均值滤波去除散在的假阳性体素，并平滑血管边缘，其中均值滤波的滤波核为3×3×3；

S2：对脑小血管区域图像进行特征提取，得到脑小血管图像特征：

采用卷积神经网络(CNN)模型(如图2所示)对滤波后的脑小血管区域图像进行特征提取的过程为：

遍历灰度值在背景上限和血管下限之间的点，取11×11×11的区域作为训练样本；

在输入层中输入训练样本，对该样本进行归一化处理，然后进行神经网络的训练；

设置第一层为卷积层，设置128个卷积核，每个卷积核大小为3×3×3，并与输入层相连进行卷积；得到128个特征map；

设置第二层为池化层，池化核大小为2×2×2，对第一层的各个map进行下采样，获得128个特征map；

设置第三层为卷积层，对第二层的各个特征map采用3×3×3大小的卷积核进行卷积，输出128个特征map；

设置第四层为Flatten层，对第三层数据进行压平；

设置第五、六、七层为全连接层对数据进行加权，其中第五层共128个类，第六层共16个类，第七层共2个类；

输出层采用softmax的激活函数，对上层卷积得到的特征进行分类，得到脑小血管图像特征；

其中CNN模型采用了2000个正样本和15000个负样本作为训练单元；并且训练集和测试集按照8:2的比例设定。

S3：基于脑小血管图像特征，对脑小血管进行追踪，获得脑小血管追踪结果；

首先采用双线性插值的方式对脑小血管区域图像进行放大，将图像中的体素(voxel)插值到各项同性的水平；在放大完成后，采用均值滤波对图像进行平滑，其中滤波核为3×3×3；

整体上采用圆柱体嵌套的方式对平滑后的放大图像中的血管进行拟合，获得拟合后的血管图像；

在本实施例中，在拟合后的血管图像上选定2个点作为初始向量，以该向量的方向为z轴进行方向计算，其中采用左手坐标系；设定初始步长为10，初始偏转角度为20°，采用随机森林方法获取圆柱体顶点位置；在训练随机森林的过程中，选取2000个正样本和46000个负样本用于随机森林的训练，在随机森林模型中，采用了5个圆柱体参数作为随机森林中的随机选取的特征，进行样本分类；

其中圆柱体参数包括：B:圆柱体内血管点数量占圆柱体体积的比值；R:圆柱体顶点在四个方向到血管边缘的距离方差；A:圆柱体中心线的偏转角度；S:圆柱体顶点周围血管点数量；H:三维Hessian矩阵的三个特征值标准差。

根据圆柱体的参数样本分类的计算结果，我们可以不断确定下一步圆柱体顶点的位置，直到无法追踪到下一步血管，此时为初步追踪结束；接下来，随机森林模型自动变化搜索角度和搜索步长，在同一步长下先按5°的倍数改变搜索角度，当搜索角度大于60°时回归初始角度并使步长+2，重复这一过程，直到当搜索半径大于基准步长的2倍仍无法找到血管时，则最终追踪结束，获得脑小血管追踪结果。

脑小血管区域还包括分支脑小血管区域；基于高连通性区域来对分支脑小血管进行搜索，所述高连通性区域是在法平面上，根据血管中心线中每一个点对应的当前圆柱体的半径变化率大于0.5的区域；其中分支脑小血管在每一步搜索到的圆柱体顶点的球形区域中进行确定，该球形区域为搜索到的圆柱体半径3倍大的球形区域；将确定的分支脑小血管添加到脑小血管追踪结果中。

S4：对脑小血管追踪结果进行筛选，获得正常脑小血管追踪结果：

首先取待筛选血管长度的1/3中每一点作为Hessian矩阵的中心点，通过Hessian矩阵的三个特征值λ₁、λ₂、λ₃判断该点是否满足管状结构；当Hessian矩阵的特征值λ₁/λ₂<20，返回λ₃的绝对值abs(λ₃)>0.05时，该点不符合管状结构。其中设定Hessian矩阵的大小为25×25×25，且三个特征值的绝对值按大小排列后λ₁＞＞λ₂，λ₃≈0。

脑小血管追踪结果除了Hessian矩阵的特征值不符合管状结构的异常脑小血管追踪结果外，还包括损失函数大于n×10，其中n是分支级数；和其他血管重复率大于30％；曲率大于80或小于5；样条曲线中30％的点在原图中灰度值小于30％的异常脑小血管追踪结果。

基于脑小血管追踪结果，建立筛选模型，将异常的脑小血管追踪结果进行筛除，获得正常的脑小血管追踪结果。

筛选结束后，使用三阶Bezier样条曲线对筛选后的脑小血管进行平滑衔接。

S5：基于正常脑小血管追踪结果，重建血管模型，输出脑小血管的结构特征。

基于平滑衔接后的正常脑小血管追踪结果，采用VTK和VMTK的算法建立三维血管模型。采用Python中VTK包的管道和VMTK包生成及渲染代码进行重建，可以使得每一支血管的量化参数会自动显示在输出的命令行窗口中。并且当模型中包括多个数据簇cluster时，如果两个cluster连接处样条曲线上的点的平均灰度值小于30％的点，则将此处的血管半径设为0，即可正确显示病灶处的情况，本实施例中给出了部分重建的三维血管模型，如图3(a)、图3(b)所示。

为了更好理解本模型，下面结合实施例对本模型做进一步地详细说明：

与目前常用的U-NET方法相比，本发明的方法拥有更高的骰子相似系数(DSC)和更低的豪斯多夫(Hausdorff)距离。特别是在小血管区域，本发明的方法显示了更好的血管分割效果。表1显示了本发明的方法与U-NET方法在DSC和Hausdorff距离上的对比。

表1

实施例中还选取了一部分的具有三级分支脑小血管区域图像作为研究图像，分别利用本发明的方法和U-NET方法对研究图像进行分割；在分割结果中，U-NET方法分割的图像只能标记血管的主干位置，无法标记血管的分支，而本发明的方法分割的图像能连贯的标记到血管的全部三级分支，由此可知，本发明的方法在目前的自动化处理当中效果是最好的，并且在血管多处中断且有非管状现象的发生(即血管狭窄导致的血流信号间断出现)，合并脑小血管对比度下降的病人当中，也可以保证二级分支被全部追踪。

在本发明的描述中，需要理解的是，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。