三维精细血管重建方法及其系统

摘要

本发明公开一种三维精细血管重建方法。所述三维精细血管重建方法包括以下步骤：步骤一、载入三维精细血管原始图像；步骤二、基于血管特性的图像预处理；步骤三、基于血管特性的特征匹配；步骤四、图像填充及边界处理。本发明同时还公开一种三维精细血管重建系统。采用本发明提供的三维精细血管重建方法及其系统，通过映射位置关系和计算整体匹配度获取修复缺失部位的填充区域，通过填充区域修复缺失部位获得完整的高精度图像，三维精细血管重建的精确性高。

说明书

技术领域

本发明涉及数字图像处理与医学成像的交叉领域，具体地，涉及一种三维精细血管重建方法及其系统。

背景技术

在生物医学领域，不同的成像技术可以观察到不同的现象和数据，可以综合考虑不同扫描设备的优势，获取真实生物软组织多种尺度下的细微观图像数据。CT图像是对人体某一部分的扫描图像，可以对血管、肿瘤等组织成像，虽然成像完整，但其缺点在于扫描数据的细节质量明显不如组织切片，而且能分辨的组织类型有限。利用光学和电子显微镜对组织切片成像，能够获取分辨清晰度更高、组织结构类型更丰富的细微观多尺度图像数据，但获得的图像数据易存在局部区域缺失的问题，因此，需要对高精度扫描图像中的局部缺失区域进行修复重建。

医学中的三维精细血管图像有着类似树形的复杂拓扑，尤其血管分支上存在分叉区域，而且分叉区域的分支数量多少不一，同时血管的粗细尺度变化大，在可视化研究中，对血管的体绘制一直是个棘手的问题，而外科手术中对主要的血管形态需要有准确的描述，以辅助医生及时选择合理治疗方法，比如从血管的三维结构上观察和确定病变位置，做出快速诊断；为手术计划提供直观的参考依据。因此一旦三维精细血管图像出现局部缺失，就会影响医生的诊断。

所以，有必要提供一种针对有缺失部位的高精度的三维精细血管图像进行重建的技术方案。

发明内容

本发明提供一种三维精细血管重建方法及其系统，通过映射位置关系和计算整体匹配度获取重建缺失部位的填充区域，通过填充区域修复缺失部位获得完整的高精度图像，三维精细血管重建的精确性高。

一种三维精细血管重建方法，包括以下步骤：

步骤一、载入三维精细血管原始图像：

所述三维精细血管原始图像包括扫描同一对象的待修复图像和参考图像，所述待修复图像为具缺失部位的血管图像，所述参考图像为完整的血管图像，所述待修复图像的精度高于所述参考图像的精度；

根据所述待修复图像指定待修复区域，并定义包含所述待修复区域的待修复块，根据所述待修复块确定在所述参考图像中位置相应、大小相同的参考映射块，所述待修复块中的待修复区域对应所述参考映射块中的参考映射区域；

步骤二、基于血管特性的图像预处理：

对载入的所述三维精细血管原始图像进行分析，确定需要增强的血管区域，对所述三维精细血管原始图像进行增强，得到三维精细血管增强图像，所述三维精细血管增强图像包括增强后的待修复图像和增强后的参考图像；

步骤三、基于血管特性的特征匹配：

根据所述三维精细血管增强图像，计算整体匹配度，得到所述参考映射块在所述待修复图像中的最优匹配块，所述参考映射块中的参考映射区域对应所述最优匹配块中的填充区域；

步骤四、图像填充及边界处理：

将所述最优匹配块的填充区域对应填充至所述待修复块的待修复区域，并对边界进行平滑处理，完成修复重建。

在本发明提供的三维精细血管重建方法的一种较佳实施例中，所述步骤二包括如下步骤：

基于Hessian矩阵对所述三维精细血管原始图像进行分析：选取Hessian矩阵对所述三维精细血管原始图像做卷积运算，求得特征值和特征向量；

建立血管区域特征函数：基于血管特性并根据所述特征值和特征向量建立血管区域特征函数；

确定需要增强的血管区域：根据所述血管区域特征函数确定需要增强的血管区域，进而对所述三维精细血管原始图像进行增强，得到三维精细血管增强图像。

在本发明提供的三维精细血管重建方法的一种较佳实施例中，所述基于Hessian矩阵对所述三维精细血管原始图像进行分析的步骤包括如下：

计算所述三维精细血管原始图像中各像素点的二阶偏导数和混合偏导数；

以计算的二阶偏导数和混合偏导数组成Hessian矩阵与所述三维精细血管原始图像做卷积运算；

求得特征值λ₁、λ₂及λ₃和与所述特征值对应的特征向量γ₁、γ₂及γ₃，且|λ₁|＜|λ₂|＜|λ₃|。

在本发明提供的三维精细血管重建方法的一种较佳实施例中，所述血管区域特征函数为：

$C\left(\lambda \right)=\left(\begin{array}{c}[1-\exp (-\frac{R_A^2}{2a^2}\left)\right][1-\exp (-\frac{R_D^2}{2d^2}\left)\right][1-\exp (-\frac{S^2}{2c^2}\left)\right],if{\lambda }_2,{\lambda }_3<0\\ 0,else\end{array}\right);$

定义R_A＝|λ₂|/|λ₃|，其中，a、c以及d分别为平面结构参数、背景区分参数以及三维血管结构参数，k₁、k₂及k₃分别为沿所述特征向量γ₁、γ₂、γ₃反方向的梯度方向的灰度变化率的大小。

在本发明提供的三维精细血管重建方法的一种较佳实施例中，所述步骤三包括如下步骤：

将所述待修复图像均分为多个与所述参考映射块大小相同的候选块，计算所述候选块和所述参考映射块的全局匹配度S_M；

计算血管段相似度S_T；

根据所述全局匹配度S_M和所述血管段相似度S_T计算整体匹配度S_ves，所述整体匹配度S_ves＝μS_M+ηS_T，其中，μ及η分别为S_M及S_T的权重系数，且μ+η＝1；

根据所述整体匹配度S_ves，得到所述参考映射块在所述待修复图像中的最优匹配块，所述参考映射块中的参考映射区域对应所述最优匹配块中的填充区域。

在本发明提供的三维精细血管重建方法的一种较佳实施例中，所述计算血管段相似度S_T的步骤包括如下步骤：

基于逐层剥取法细化步骤二中所述需要增强的血管区域，提取单像素中心线，由所述三维精细血管增强图像得到三维精细血管二值图，所述三维精细血管二值图包括所述待修复图像的二值图和所述参考图像的二值图；

根据所述待修复图像的二值图提取所述候选块的端点集，根据所述参考图像的二值图提取所述参考映射块的端点集，分别求取端点特征向量；

根据所述端点特征向量计算血管段相似度S_T。

在本发明提供的三维精细血管重建方法的一种较佳实施例中，所述步骤四包括如下步骤：

将所述最优匹配块的填充区域对应填充至所述待修复块的待修复区域；

将所述填充区域的边界进行插值迭代，不断逼近直至收敛，将所述填充区域、以及所述待修复区域的边界平滑过渡，完成修复重建。

本发明还提供一种三维精细血管重建系统，包括图像载入模块，用于载入三维精细血管原始图像，所述三维精细血管原始图像包括扫描同一对象的待修复图像和参考图像，所述待修复图像为具缺失的血管图像，所述参考图像为完整的血管图像，所述待修复图像的精度高于所述参考图像；并根据所述待修复图像指定待修复区域，并定义包含所述待修复区域的待修复块，根据所述待修复块确定在所述参考图像中位置相应、大小相同的参考映射块，所述待修复块中的待修复区域对应所述参考映射块中的参考映射区域；

图像预处理模块，用于对所述图像载入模块载入的所述三维精细血管原始图像进行分析，确定需要增强的血管区域，对所述三维精细血管原始图像进行增强，得到三维精细血管增强图像，所述三维精细血管增强图像包括增强后的待修复图像和增强后的参考图像；

特征匹配模块，用于根据所述图像预处理模块得到的三维精细血管增强图像，计算整体匹配度，得到所述参考映射块在所述待修复图像中的最优匹配块，所述参考映射块中的参考映射区域对应所述最优匹配块中的填充区域；

填充处理模块，用于将所述特征匹配模块得到的所述最优匹配块的填充区域对应填充至所述待修复块的待修复区域，并对边界进行平滑处理，完成修复重建。

在本发明提供的三维精细血管重建系统的一种较佳实施例中，所述图像预处理模块包括：

第一运算单元，选取Hessian矩阵对所述三维精细血管原始图像做卷积运算，求得特征值和特征向量；

函数创建单元，根据所述第一运算单元求得的特征值和特征向量并基于血管特性建立血管区域特征函数；

血管增强单元，根据所述函数创建单元建立的血管区域特征函数确定需要增强的血管区域，进而对所述三维精细血管原始图像进行增强，得到三维精细血管增强图像。

在本发明提供的三维精细血管重建系统的一种较佳实施例中，所述特征匹配模块包括第一计算单元、第二计算单元及第三计算单元，分别用于计算全局匹配度、血管段相似度及整体匹配度。

相较于现有技术，本发明提供的三维精细血管重建方法及其系统，具有以下有益效果：

一、本发明提供的三维精细血管重建方法对需要增强的血管区域进行了特征值提取，将三维精细血管原始图像进行了图像增强，便于对源于同一血管样本的待修复图像和参考图像进行结构匹配，提高了三维精细血管重建的效率，并能够保证精度；

二、本发明重新定义了基于Hessian矩阵的血管区域特征函数，考虑了三维精细血管原始图像的灰度信息，可以有效提高图像增强的效果；同时考虑了血管的分叉区域以及弧度较大的区域等个体差异性的影响，能够很好的解决个体差异性的问题，即使在特殊情况也不会判断错误，具有更广泛的适用性；

三、本发明提供的三维精细血管重建方法基于血管特性而成，考虑到血管结构的管状特性和较强的自相似性，从而指导源于同一血管样本的待修复图像和参考图像的匹配，提高了血管结构特征匹配的可靠性和准确性，进一步增大三维精细血管重建的精度；

四、本发明通过待修复图像和参考图像的匹配，得到参考映射块在待修复图像中的最优匹配块，对应填充之后再对边界进行平滑处理，修复精度高。

附图说明

图1为本发明提供的三维精细血管重建方法的修复原理图；

图2为本发明提供的三维精细血管重建方法的流程图；

图3(a)为本发明提供的三维精细血管原始图像中的参考图像；

图3(b)为本发明提供的三维精细血管原始图像中的待修复图像；

图4为图2所示三维精细血管重建方法中步骤S2的流程图；

图5为图4所示步骤S2中步骤S21的流程图；

图6(a)为图3(a)所示增强后的参考图像；

图6(b)为图3(b)所示增强后的待修复图像；

图7为图2所示三维精细血管重建方法中步骤S3的流程图；

图8为图7所示步骤S3中步骤S32的流程图；

图9(a)为图6(a)所示细化后的参考图像；

图9(b)为图6(b)所示细化后的待修复图像；

图10为图2所示三维精细血管重建方法中步骤S4的流程图；

图11(a)为待修复的一种三维精细血管结构图；

图11(b)为图11(a)所示修复后的三维精细血管结构图；

图12为本发明提供的三维精细血管重建系统的结构框图；

图13为图12所示图像预处理模块的结构框图；

图14为图12所示特征匹配模块的结构框图；

图15为图12所示填充处理模块的结构框图。

具体实施方式

下面将结合附图和实施方式对本实用新型(发明)作进一步说明。请参阅图1，为本发明提供的三维精细血管重建方法的修复原理图。本发明需解决的技术问题是修复具缺失部位的三维精细血管图像，记为待修复图像1。本发明还提供了用于修复所述待修复图像1的参考图像3。具体地，在不破坏血管样本的前提下，为了获得血管样本的微观多尺度图像数据，对一份血管样本进行序列连续切片，采用分辨率高的光学和电子显微镜进行扫描获得所述待修复图像1，由于分辨率越高，能扫描样本区域的范围越小，对血管样本进行细微观多尺度扫描获得的图像数据存在缺失部位，即所述待修复图像1为具缺失部位的高精度血管图像；采用显微CT/MRI和同步辐射CT直接扫描同一血管样本，获得所述参考图像3，所述参考图像3为完整的低精度血管图像；此低精度和高精度是相对概念，即所述待修复图像1的扫描精度高于所述参考图像3的扫描精度。

T1、指定所述待修复图像1中所述缺失部位为待修复区域10，定义从所述待修复区域10向外延伸获得的立方体扩张区域为待修复块12；

T2、通过映射位置关系在所述参考图像3中获得与所述待修复块12位置相应、大小相同的参考映射块32，所述待修复块12中的待修复区域10对应所述参考映射块32中的参考映射区域30；

T3、在所述待修复图像1中搜索匹配与所述参考映射块32相似的最优匹配块52，再通过映射位置关系获得与所述待修复区域10相对应的填充区域50；

T4、将所述填充区域50对应填充至所述待修复区域10，所述填充区域50修复所述缺失部位，对所述填充区域50和所述待修复区域10的边界进行平滑处理，得到自然过渡的纹理图像，完成修复重建，得到完整的高精度三维精细血管图像。

请参阅图2，为本发明提供的三维精细血管重建方法的流程图。本发明提供一种三维精细血管重建方法，包括如下步骤：

步骤S1、载入三维精细血管原始图像：

所述三维精细血管原始图像包括扫描同一对象的待修复图像1和参考图像3，所述待修复图像1为具缺失部位的血管图像，所述参考图像3为完整的血管图像，所述待修复图像1的精度高于所述参考图像3的精度。

在本实施例中，具体地，请参见图3(a)及图3(b)，其中，图3(a)为本发明提供的三维精细血管原始图像中的参考图像；图3(b)为本发明提供的三维精细血管原始图像中的待修复图像。需要说明的是，为了更清楚的描述本发明提供的方法中的关键步骤，本发明对上述源于同一血管样本的参考图像和待修复图像进行了适当截取，并以此为实例进行说明，但不局限于上述公开的图像。

根据所述待修复图像1指定待修复区域10，并定义包含所述待修复区域10的待修复块12，根据所述待修复块12确定在所述参考图像3中位置相应、大小相同的参考映射块32，所述待修复块12中的待修复区域10对应所述参考映射块32中的参考映射区域30。

在本实施例中，指定所述待修复图像1中所述缺失部位为待修复区域10，定义从所述待修复区域10向外延伸获得的立方体扩张区域为待修复块12。通过映射位置关系在所述参考图像3中获得与所述待修复块12位置相应、大小相同的参考映射块32，所述待修复块12中的待修复区域10对应所述参考映射块32中的参考映射区域30。

步骤S2、基于血管特性的图像预处理：

对载入的所述三维精细血管原始图像进行分析，确定需要增强的血管区域，对所述三维精细血管原始图像进行增强，得到三维精细血管增强图像，所述三维精细血管增强图像包括增强后的待修复图像1和增强后的参考图像3。

步骤S3、基于血管特性的特征匹配：

根据所述三维精细血管增强图像，计算整体匹配度，得到所述参考映射块32在所述待修复图像1中的最优匹配块52，所述参考映射块32中的参考映射区域30对应所述最优匹配块52中的填充区域50。

步骤S4、图像填充及边界处理：

将所述最优匹配块52的填充区域50对应填充至所述待修复块12的待修复区域10，并对边界进行平滑处理，完成修复重建。

请参阅图4，为图2所示三维精细血管重建方法中步骤S2的流程图。由于采集设备光强的误差，三维精细血管等医学图像中感兴趣部位(即血管区域)的灰度值可能并不明显，为了便于对图像进行修复，需要先对图像进行增强处理，增强关键结构的匹配。

根据结构张量能够识别图像中的边缘、角点以及平坦区域，结构张量是对图像中的像素利用矩阵组织的数据结构，其形式便是Hessian矩阵。对于形状特征特殊的血管图像，血管截面的灰度成像高斯分布的特性，可以利用结构张量有效得到血管所需特征，进而对其进行特征增强。

所述基于血管特性的图像预处理步骤S2包括：

步骤S21、基于Hessian矩阵对所述三维精细血管原始图像进行分析：选取Hessian矩阵对所述三维精细血管原始图像做卷积运算，求得特征值和特征向量。

请同时参阅图5，为图4所示步骤S2中步骤S21的流程图。

所述基于Hessian矩阵对所述三维精细血管原始图像进行分析的步骤S21包括：

步骤S21-1、计算所述三维精细血管原始图像中各像素点的二阶偏导数和混合偏导数；

步骤S21-2、以计算的二阶偏导数和混合偏导数组成Hessian矩阵与所述三维精细血管原始图像做卷积运算；

步骤S21-3、求得特征值λ₁、λ₂及λ₃和与所述特征值对应的特征向量γ₁、γ₂及γ₃，且|λ₁|＜|λ₂|＜|λ₃|。

具体过程为：

设V_x,V_y,V_z分别为x,y,z方向的偏导数，计算所述三维精细血管原始图像中各像素点的二阶偏导数和混合偏导数，其中V_σ＝G_σ*V，*代表卷积算子，通过高斯函数(1)获得不同尺度σ中的线性增强滤波。进行局部特性分析时，以当前处理像素点为中心，在当前所处理的图像数据上，取半宽为3σ的矩形窗口，足以包含血管的直径。

所述高斯函数(1)为：

结构张量(是一个三维的Hessian矩阵)定义如下式(2)：

其中是笛卡尔内积；存在一个三维正交矩阵S，使得其中Λ＝diag(λ_m)的对角元素是的特征值，组成矩阵S的各行向量即为结构张量Υ的特征向量，分别记为λ₁、λ₂及λ₃，且满足|λ₁|＜|λ₂|＜|λ₃|，与三个特征值分别一一对应的三个特征向量记为γ₁、γ₂及γ₃，其中γ₁⊥λ₂,γ₁⊥λ₃,γ₂⊥λ₃。结构张量Υ的特征值λ₁、λ₂及λ₃能够体现出图像灰度在相应特征向量方向上的变化速度。图像的边缘强度大，可以量化为结构张量Υ较大的特征值。

根据三个特征值|λ₁|＜|λ₂|＜|λ₃|，对矩阵的特征值及对应的形状结构进行如下分析：

a.|λ₁|≈0,|λ₂|≈0,|λ₃|≈0,可判定为平坦区域；可以判定为面状结构；可以判定为管状结构；可判定为球状结构。

通过上述分析，对于理想的三维精细血管图像，其特征值及特征向量的关系应该如下式(3)：

其中λ₂,λ₃均为负数，表示了在三维精细血管图像中感兴趣部位为高亮区而背景为阴暗区；矩阵的两个特征向量γ₂,γ₃组成的面所代表的是三维精细血管的切面，而最小特征值所对应的特征向量γ₁代表的是三维精细血管的走向，即三维精细血管的延伸方向。

步骤S22、建立血管区域特征函数：基于血管特性并根据所述特征值和特征向量建立血管区域特征函数。

现有的血管区域增强算法：根据所述特征值λ₁、λ₂、λ₃和所述特征向量γ₁、γ₂、γ₃，首先定义了变量R_A、R_B及S，分别如下式(4)、(5)、(6)：

R_A＝|λ₂|/|λ₃|(4)；

变量R_A的值可以估计当前处理图像属于管状结构的可能性，其值越大，越可能属于管状结构；

$R_B=\left|{\lambda }_1\right|/\sqrt{{\lambda }_2{\lambda }_3}---\left(5\right);$

变量R_B的值可以估计当前处理图像属于盘状结构的可能性，其值越大，越可能属于盘状结构；

$S=\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}---\left(6\right);$

变量S取所述三个特征值的平方和，可以模糊判断当前处理图像属于结构性较强的区域。

根据所述变量R_A、R_B及S定义的现有血管区域函数为下式(7)：

$C\left(\lambda \right)=\left(\begin{array}{c}[1-\exp (-\frac{R_A^2}{2a^2}\left)\right]\left[\exp \right(-\frac{R_B^2}{2b^2}\left)\right][1-\exp (-\frac{S^2}{2c^2}\left)\right],if{\lambda }_2,{\lambda }_3<0\\ 0,else\end{array}\right)---\left(7\right);$

其中，a、b以及c分别为平面结构参数、球状结构参数及背景区分参数，根据经验以及所选尺度选取。

针对真实的三维精细血管图像，血管粗细及拓扑都不规则，血管粗细不等，对于非常细的血管，其交叉点几乎可以忽略，但是对于较粗的血管，其分叉区域就不容忽视了。

根据公式(5)计算变量R_B时，若当前处理图像为盘状结构或分叉区域，λ₁是和λ₂、λ₃大小相当的负数，这两种情况下，变量R_B的值均较大；若当前处理图像为血管弧度较大的区域，R_B的值也会较大。

但根据公式(7)计算现有血管区域函数时，R_B越小，当前处理图像属于血管区域的权值越大，因而在特殊情况下，根据公式(7)判断当前处理图像会出现错误，所以应该找回被忽略的分叉区域和血管弧度较大的区域。

本发明假设所述特征向量γ₁、γ₂、γ₃的反方向分别对应为τ₁,τ₂,τ₃，所述特征向量γ₁、γ₂、γ₃指向血管壁，在一定概率下，τ₁,τ₂,τ₃中有至少一个方向指向血管延伸方向，如果沿着τ₁,τ₂,τ₃其中某一方向像素灰度值变化缓慢，则可以判定为有效的血管区域。

定义沿方向τ的梯度方向对于尺度σ，计算得到沿此梯度方向的变化率的大小k，设三个方向τ₁,τ₂,τ₃的灰度变化率大小分别为k₁,k₂,k₃，本发明定义新变量即血管角点区域特性R_D，如下式(8)：

$R_D=1-k_1{k}_2{k}_3/\sqrt{{\left({\lambda }_2{\lambda }_3\right)}^3}---\left(8\right).$

根据所述变量R_A、R_B及新变量R_D，定义新的血管区域函数为下式(9)：

$C\left(\lambda \right)=\left(\begin{array}{c}[1-\exp (-\frac{R_A^2}{2a^2}\left)\right][1-\exp (-\frac{R_D^2}{2d^2}\left)\right][1-\exp (-\frac{S^2}{2c^2}\left)\right],if{\lambda }_2,{\lambda }_3<0\\ 0,else\end{array}\right)---\left(9\right);$

其中，a、c以及d分别为平面结构参数、背景区分参数以及三维血管结构参数，根据经验以及所选尺度选取。

若当前处理图像为分叉区域，k₁、k₂及k₃中只要任意一个值较小，则变量R_D的值较大；若当前处理图像为血管弧度较大的区域，k₂及k₃与λ₂及λ₃大小相当，k₁的值相对较小，则变量R_D的值较大；若当前处理图像为典型的直血管区域，k₁≈0，变量R_D的值依然较大；因此新变量R_D适用多种血管区域，不会对其他特征项产生干扰。

步骤S23、确定需要增强的血管区域：根据所述血管区域特征函数确定需要增强的血管区域，进而对所述三维精细血管原始图像进行增强，得到三维精细血管增强图像。所述三维精细血管增强图像包括增强后的待修复图像1和增强后的参考图像3。

在本实施例中，具体地，请参见图6(a)及图6(b)，其中，图6(a)为图3(a)所示增强后的参考图像；图6(b)为图3(b)所示增强后的待修复图像。对比发现利用本发明提供的血管区域特征函数进行图像增强具有较好的效果。

请参阅图7，为图2所示三维精细血管重建方法中步骤S3的流程图。所述基于血管特性的特征匹配步骤S3包括如下步骤：

步骤S31、将所述待修复图像1均分为多个与所述参考映射块32大小相同的候选块，计算全局匹配度S_M。

根据所述三维精细血管增强图像，即所述待修复图像1和所述参考图像3，对两幅图像中血管的整体分布进行相似度匹配，计算所述候选块和所述参考映射块32的全局匹配度S_M。

新建0-1矩阵B₁，将所述待修复图像1的所述候选块中属于血管的相应位置设为1，属于背景的相应位置设为0；新建0-1矩阵B₂，将所述参考图像3的所述参考映射块32中属于血管的相应位置设为1，属于背景的相应位置设为0；对所述矩阵B₁和所述矩阵B₂相应位置上进行“异或运算”，得到结果矩阵M_R，所述矩阵M_R中“0”越多，即相同之处越多，则说明所述待修复图像1的所述候选块和所述参考图像3的所述参考映射块32从轮廓上看越相似。

定义全局匹配度S_M为下式(10)：

$S_M=\frac{1}{n\times n\times n}\underset{(x,y,z)\in B_1,B_2}{\Sigma }u(x,y,z)---\left(10\right);$

其中，定义u(x,y,z)为下式(11)：