颅骨组织工程Mosaic结构血管网参数化设计方法

摘要

本发明公开了一种颅骨组织工程Mosaic结构血管网参数化设计方法，包括步骤一、板障静脉血管特征提取，颅骨模型重建后，通过模型预处理，比较空腔结构半径大小，区分出多孔结构与静脉血管结构，然后，提取板障静脉血管特征，所述板障静脉血管特征包括拓扑特征和几何特征，所述拓扑特征反映血管在空间中的分布，所述几何特征包括血管分支角度和血管半径；步骤二、血管结构参数化模型设计，首先，Mosaic结构参数化设计，然后，定义血管路径；最后生成具有多级特征的参数化血管模型。本发明中，通过调整Mosaic结构参数，可以调整血管网疏密程度，并可针对某一区域进行血管结构调整。能使血管均匀、充分地分布在设定空间中，作为颅骨血管支架是可行的。

说明书

技术领域

本发明属于组织工程技术领域，本发明涉及一种颅骨组织工程Mosaic结构(马赛克结构)血管网参数化设计方法。

背景技术

骨组织工程技术研制出的具有生命活性的人工骨替代物，可对损伤的颅骨组织进行修复和重建，从而维持或改善颅骨的功能和结构，在临床上具有很好的应用前景。但在颅骨多孔支架的体外细胞整合研究中，发现骨细胞的存活率低，成骨缓慢，特别是对大面积支架，在植入的第1周中心处就有大量的细胞死亡，其主要原因是未在支架中考虑血管化问题，支架中没有血管结构，细胞的供血和供氧受到限制，以致骨细胞不能深入支架内部，造成成骨能力下降。

解剖学研究表明，在颅骨外板和内板中间的板障空间中存在有血管组织，称为板障静脉。因此，为了给骨细胞生存提供充足的营养物质和氧气，理想的颅骨支架应具有与天然板障静脉相似的血管结构，才能实现血液与代谢物的传输，从而提高细胞存活率。对于血管结构，现有研究大多数是针对人体软组织部位，例如肺部血管、肝脏血管。而颅骨板障静脉，由于它所处的位置特殊，形态复杂，尺度微小，目前还没有研究人员进行该组织血管的设计工作。

具有血管结构的颅骨支架对细胞的生长很重要，为此本发明通过显微CT扫描天然颅骨样本，重建板障中的静脉血管结构，提取静脉血管空间结构、分叉、管径等形态特征，据此建立血管的参数化模型，可快速、灵活地对血管的结构进行设计，使其适用于任意大小的支架，可为颅骨组织工程支架血管化的实现奠定基础。

发明内容

针对颅骨组织工程支架由于缺乏血管结构从而限制了细胞的生长及成骨量的形成等问题，提出了一种基于人体颅骨板障静脉血管特征的Mosaic结构血管网参数化设计方法。

本发明提出的一种颅骨组织工程Mosaic结构血管网参数化设计方法，包括以下步骤：

步骤一、板障静脉血管特征提取，包括：

1-1)颅骨模型重建，对一颅骨样本进行显微CT扫描；采用三维可视化软件Avizo对扫描数据进行颅骨模型三维重建；对该颅骨模型采用阈值分割，提取板障中的空腔结构，该空腔结构至少包括骨组织间的孔洞和静脉血管；

1-2)模型预处理，比较空腔结构半径大小，区分出多孔结构与静脉血管结构；

采用Avizo软件提取所述空腔结构的中心线，进而计算统计出自该中心线到空腔结构边界的距离r_i即半径，在所述空腔结构中，将r_i＜0.4mm的结构定义为多孔结构，将r_i≥0.4mm的结构为静脉血管结构；按照静脉血管结构半径r_i的不同将血管结构分为1～4级，自1级血管结构至4级血管结构依次为高级血管结构至低级血管结构；所述静脉血管结构对应的中心线为板障静脉血管中心线；

1-3)提取板障静脉血管特征，所述板障静脉血管特征包括拓扑特征和几何特征，所述拓扑特征反映血管在空间中的分布，所述几何特征包括血管分支角度和血管半径；其中：

提取拓扑特征的过程是：将板障静脉血管中心线投影到血管网具有最大投影面积的平面，进行曲线直线化处理后，在该投影中，低级别的静脉血管结构镶嵌在高级别的静脉结构中，且所述静脉血管结构在投影区域内均为封闭多边形，将该封闭多边形定义为Mosaic结构；Mosaic结构的面积范围为0.3～1256mm²；

提取几何特征的过程是：通过抽样统计得出血管分支角度α＝80°～100°；血管横截面内接圆半径r_i的范围为0.4～1.3mm，所述血管横截面内接圆半径r_i即为血管半径；

步骤二、血管结构参数化模型设计，包括：

2-1)Mosaic结构参数化设计，设：颅骨空间大小为X＝a×b×c，依据步骤1-3)中板障静脉血管特征，在平面区域大小为a×b内利用power图生成第一级Mosaic结构参数化模型；其中power图中控制圆的面积大小和比例与步骤1-3)中Mosaic结构的面积大小和所占比例一致；在第一级Mosaic结构参数化模型的每个Mosaic结构中再利用power图生成第二级Mosaic结构，直至生成最高4级相互嵌套的Mosaic结构；

对上述Mosaic结构参数化模型进行干涉判断，即，当两节点间距离小于1.2倍血管半径时，则认为生成的新节点为干涉点，将其删除；

2-2)血管路径的定义，内容如下：

①三维Mosaic结构的生成：

采用分层逼近法将步骤2-1)获得的Mosaic结构参数化模型转换为三维空间结构，即为三维Mosaic结构；过程是：以Mosaic结构参数化模型所在平面P的中心点O_c为圆心，在所述平面P的垂直方向做n_c个同心圆，设定相邻同心圆之间的半径差为z_c；根据颅骨样本空间弯曲曲率R_c，平面区域大小a×b，则半径差z_c为：

对位于半径最小的同心圆内的节点，在范围[0,z_c]内随机赋予Z坐标；自半径最小至半径最大的同心圆内的节点，分别在范围[z_c,2×z_c]，[z_c,3×z_c]，……，[z_c,n_c×z_c]中随机赋予Z坐标；同时保证X,Y坐标相同的节点，Z坐标也保持相同，从而完成了Mosaic结构从二维到三维的转换；

②血管路径生成：

将上述得到的三维Mosaic结构的顶点定义为血管节点，将拟定的血流流向定义为血流矢量v_ci，将处于边界的血管节点作为血管起始点；定义当前节点指向下一节点的矢量为v_a-b，依据血流矢量v_ci与矢量v_a-b的最小夹角确定下一节点，依次生成血管路径；

③血管中心线参数化模型建立：

通过B样条曲线对上述生成的血管路径进行拟合，定义拟合曲线中第一点为起始点，最后一点为末端点，若拟合曲线起始点位于高一级血管路径上，则保证拟合曲线起始点切矢与高一级血管路径夹角在80°～100°之间，且除端点之外的其它节点均满足二阶连续可导，即得血管中心线参数化模型；

2-3)参数化血管模型生成，

结合上述建立的血管中心线参数化模型以及步骤1-3)中提取的几何特征，并采用三维制图软件对各级血管进行批量扫掠，完成各级血管求和布尔运算后根据设定的血管壁厚进行抽壳操作，即可生成具有多级特征的参数化血管模型。

进一步讲，在步骤1-2)模型预处理中，按照静脉血管结构半径r_i的不同将血管结构分为如下4级：

本发明中，所述颅骨空间尺寸X＝a×b×c＝100mm×100mm×3mm，空间曲率半径R_c为125mm。

血管壁厚与血管半径之比为相对厚度，取值范围为0.04～0.1。

与现有技术相比，本发明的有益效果是：

为了验证本发明方法的可行性，以血管空间尺寸X＝100mm×100mm×3mm为例，构建两组不同参数的血管实例，来评价本发明的效果。

实例一：具有四级血管结构的血管网参数取值如表1所示，设计结果如图1(a)和图1(b)所示。该模型中主血管体积为1074mm³，血管总体积为2877mm³，该空间内血管密度为1.6％。

表1血管设计参数

实例二：为了表明Mosaic结构面积S大小对血管网形态的影响，在合理范围内将S减小至8mm²-1700mm²，其他参数如实例一保持不变，设计结果如图1(c)和图1(d)所示。实例二中血管总体积为4270mm³，空间内血管密度为2.1％。

由上述两实例可看出，两组参数所得到的血管网分布均匀，小血管也都能充分地填充主血管之间的空隙，有效地提高了空间中的血管密度。通过调整Mosaic结构参数，可以调整血管网疏密程度，并可针对某一区域进行血管结构调整。由此可见，该Mosaic结构血管网参数化设计方法能使血管均匀、充分地分布在设定空间中，作为颅骨血管支架是可行的。

附图说明

图1(a)是血管模型设计实例一；

图1(b)是血管模型设计实例一侧视图；

图1(c)是是血管模型设计实例二；

图1(d)血管模型设计实例二侧视图；

图2是血管设计流程图；

图3是颅骨样本；

图4(a)是颅骨扫描切片图；

图4(b)是颅骨样本三维模型；

图5是板障中的空腔结构；

图6(a)是板障中空腔结构中心线；

图6(b)是板障静脉血管中心线投影图；

图6(c)是四级血管镶嵌示意图；

图6(d)是血管局部特征示意图；

图7是二维平面power图；

图8是具有二级Mosaic结构的平面血管模型；

图9(a)是Mosaic结构坐标转换示意图；

图9(b)是转换后的三维Mosaic结构；

图10是血管路径生成示意图；

图11(a)是血管中心线模型；

图11(b)是参数化血管模型；

图11(c)是参数化血管模型侧视图。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。

本发明提出的一种颅骨组织工程Mosaic结构血管网参数化设计方法，设计流程如图2所示，包括以下步骤：

步骤一、板障静脉血管特征提取，包括：

1-1)颅骨模型重建，选取人体自然颅骨额骨部位(由天津医科大学提供)，尺寸约100mm×70mm×3mm，如图3所示。通过显微CT扫描(天津三英精密仪器有限公司所提供的X射线三维显微镜nanoVoxel设备)，得到图像共计720帧，其图像分辨率为49.8μm，扫描切片图如图4(a)所示。运用三维可视化软件Avizo对扫描数据进行颅骨模型三维重建，获得颅骨样本三维模型如图4(b)所示，对重建的体数据进行阈值分割，提取板障中的空腔结构，结果如图5所示。所述空腔结构至少包括骨组织间的孔洞和静脉血管。

1-2)模型预处理，比较1-1)所述空腔结构半径大小，区分出多孔结构与静脉血管结构。预处理过程如下：①采用Avizo软件提取所述空腔结构中心线，如图6(a)所示，进而计算统计出自该中心线到空腔结构边界的距离r_i即半径；将r_i＜0.4mm的结构定义为多孔结构，将r_i≥0.4mm的结构定义为静脉血管结构；统计不同半径范围在整个空腔结构中所占的比例，统计结果如表2所示。按照静脉血管结构半径r_i的不同将血管结构分为1～4级，自1级血管结构至4级血管结构依次为高级血管结构至低级血管结构，如表3所示；所述静脉血管结构对应的中心线为板障静脉血管中心线；血管半径的统计值可为血管设计提供参数依据。

表2空腔结构半径范围统计结果

表3

1-3)提取板障静脉血管特征，所述板障静脉血管特征包括拓扑特征和几何特征，所述拓扑特征反映血管网在空间中的分布，所述几何特征反映血管分支角度和血管半径；其中：

提取拓扑特征的过程是：将板障静脉血管中心线投影到血管网具有最大投影面积的平面，进行曲线直线化处理，在该投影中，低级别的静脉血管结构镶嵌在高级别的静脉结构中，如图6(b)所示。定义静脉血管结构是由多级不同直径的血管层层叠加而成，如图6(c)所示。且所述静脉血管结构在投影区域内均为封闭多边形，将该封闭多边形定义为Mosaic结构；通过对Mosaic结构的面积进行统计，其面积范围为0.3～1256mm²；，面积比例如表4所示。

表4颅骨样本中Mosaic结构面积及比例

提取几何特征的过程是：通过抽样统计得出血管分支角度α＝80°～100°，血管横截面内接圆半径r_i的范围为0.4～1.3mm，所述血管横截面内接圆半径r_i即为血管半径，如图6(d)所示。

步骤二、血管结构参数化模型设计，设：

颅骨空间尺寸X＝a×b×c＝100mm×100mm×3mm，空间曲率半径R_c为125mm，血管结构参数化模型设计过程包括：

2-1)Mosaic结构参数化设计，

静脉血管结构是由不同级别的Mosaic结构镶嵌而成，需要对Mosaic结构面积进行参数化定义；本发明采用power图来生成Mosaic结构参数化模型。

设S＝{p₁,p₁,...,p_n}为二维欧式空间平面上的点集，对每个生成元p_i附以权重w_i(i＝1,2,…,n)。用power距离

d_w(p,p_i)＝||p-p_i||-w_i>

定义的Vor(p,d_pw)＝Vor_pw＝{V(p_i),...,V(p_n)}称为由集合S生成的加权power图，或简称为S的power图，称集合V(p_i)为与p_i相关联的power多边形。图7所示为二维平面上的power图，定义圆点为生成元p_i，线条为power边，生成元p_i生成的圆为控制圆(用来控制Mosaic结构面积大小)。Mosaic结构参数化模型生成过程如下：

设：颅骨空间大小为X＝a×b×c，颅骨样本空间弯曲曲率为R_c；在平面a×b中生成生成第一级Mosaic结构参数化模型；在血管设计过程中，power图中控制圆的大小和比例与步骤1-3)中Mosaic结构的面积大小、所占比例一致；

将颅骨样本板障静脉中心线投影到血管网具有最大投影面积的平面，通过计算机视觉识别计算Mosaic结构面积，Mosaic结构面积S和所占比例t_p如表4所示。令控制圆面积s_c等于Mosaic结构面积S，则控制圆个数c_p可由公式(2)确定。

通过表4与公式(2)中的Mosaic结构面积、比例以及平面大小a×b，可以确定控制圆面积大小与个数，生成相应的power图，即Mosaic结构。以二级Mosaic结构为例进行说明：首先依据公式2生成第一级控制圆和power图，然后在其每一个Mosaic结构中依据表4统计数据生成下一级控制圆和power图，从而生成具有二级Mosaic结构的平面血管模型，如图8所示。需要说明，该方法可生成具有多级Mosaic结构的平面血管模型，但依据静脉血管在颅骨中的分布密度，建议Mosaic结构分级最多分为4级。生成Mosaic结构的平面血管模型后，需要对上述Mosaic结构参数化模型进行干涉判断，即，当两节点间距离小于1.2倍血管半径时，则删除其中一个节点。

2-2)血管路径的定义，内容如下：

采用power图生成Mosaic结构后，需要将二维Mosaic结构转换为三维Mosaic结构，通过扫掠、抽壳操作生成血管结构。

①三维Mosaic结构的生成

采用分层逼近法将步骤2-1)所述Mosaic结构转换为三维Mosaic结构，方法如下：找出平面a×b中心点O_c，以O_c为圆心，在Z轴方向(平面垂直方向)做n_c个同心圆，定义相邻同心圆之间的半径差为z_c，如图9(a)所示，其中n_c值越大，血管空间弯曲越平滑，弯曲效果如图9(b)所示。根据颅骨样本空间弯曲曲率R_c，最大偏置距离z_c可通过公式3确定。假定设计的平面区域大小为a×b，则最大偏置距离

位于半径最小的同心圆内的节点，在范围[0,z_c]内随机赋予Z坐标，自半径最小至半径最大的同心圆内的节点，分别在范围[z_c,2×z_c]，[z_c,3×z_c]，……，[z_c,n_c×z_c]中随机赋予Z坐标；同时保证X,Y坐标相同的节点，Z坐标也保持相同，从而完成了Mosaic结构从二维到三维的转换。

②血管路径生成：

当前的三维Mosaic结构中需要确定血管生成路径，并保证血管分支角度满足统计要求，便于后续扫掠操作；以二级血管为例说明路径生成过程；在图10中，粗线条和细线条分别表示第一、二级血管中心线，粗箭头和细箭头代表第一、二级血流矢量，即血流方向，(N_b1…N_b10)、(N_s1…N_s10)分别为第一、二级血管中的节点。

首先，生成第一级血管路径。筛选出第一级血管节点N_bi，如图10所示，设定节点N_b1为起始节点；设定第一级血流矢量v_c1竖直向上，矢量v_c1用来筛选血管路径中的节点，同时表示第一级血管中整体的血流方向。定义当前节点指向下一节点的矢量为v_a-b，(以图10中节点N_b2为例，N_b9与N_b3均与N_b2相连,则由N_b2指向N_b9与N_b3的矢量分别为v_b2-b9与v_b2-b3)矢量v_c1与矢量v_a-b夹角最小(如图10中v_b2-b9与v_c1夹角最小，为52°)的节点(如图10中N_b9)作为下一节点，依次类推，直到第一级血管所有节点分配完毕。

其次，生成第二级血管路径。筛选出第二级血管中节点N_si，如图10所示。设定位于第一级血管中的节点(如图10中N_s5)为起始节点，定义垂直于主血管的矢量v_c2为第二级血流矢量，定义当前节点指向下一节点的矢量为v_a-b，与第一级血管路径生成方法类似，矢量v_c2与矢量v_a-b夹角最小的节点作为下一节点，直至所有节点分配完毕。表5显示了血管路径生成过程。

表5血管路径生成过程

③血管中心线参数化模型建立：

当前血管路径为空间折线，可通过B样条曲线进行拟合，以便完成扫掠操作。设：血管路径中共有n个节点，定义拟合曲线中第一点为起始点p₀，最后一点为末端点p_n-1，V_s为血管路径起始点P₀的切矢，V_e为血管路径末端节点P_n-1的切矢，V_m为主血管中的血流方向。如果起始点P₀位于高一级血管上，则设定V_s与V_m呈α角(参考步骤1-3)，α夹角范围为80°～100°)，否则起始切矢V_s由血管路径起始两节点(P₀,P₁)确定；末端切矢V_e由最末端两节点确定(P_n-2,P_n-1)。结合初始切矢V_s、末端切矢V_e以及血管路径节点，采用B样条曲线完成血管路径拟合，每条曲线在除端点之外的其它节点均满足二阶连续可导，即可获得血管中心线模型，如图11(a)所示。

2-3)曲面扫掠参数化定义

结合血管中心线参数化模型以及步骤1-3)中提取的几何特征，并采用三维制图软件对各级血管进行批量扫掠，各级血管半径可由表2中统计数据得到，在上述建立的血管中心线参数化模型起始点出生成扫掠截面，沿着中心线模型进行扫掠，即可得到不同半径的血管实体模型，完成各级血管实体模型求和布尔运算后，根据设定的血管壁厚(本发明中，血管壁厚与血管半径之比为相对厚度，取值范围为0.04～0.1)进行抽壳操作，即可生成具有多级特征的参数化血管模型，如图11(b)和图11(c)所示。

2.最佳实施例：

(1)样本数据获取

本发明通过采用显微CT获取颅骨板障静脉，并通过Avizo软件进行阈值调节，重建颅骨模型，通过AutoSkeleton命令提取血管中心线，统计血管半径信息，将重建模型导出为STL格式。对提取血管中心线投影到平面，基于Visual Studio平台采用OpenCV库编译图像识别程序，对其中封闭区域像素点个数进行计算，获取封闭区域面积。通过Geomagic软件对STL重建模型进行局部特征提取，获取分支角度，截面特征。

(2)2D Mosaic拓扑结构生成

本发明采用Matlab程序生成2D Mosaic拓扑结构，已完成相关代码编写。依据统计的Mosaic结构面积，逐级生成各级控制圆与Power图，并将Power图导出为DXF格式。

(3)3D中心线模型生成

本发明采用Matlab用分层逼近法实现2D Mosaic结构向3D Mosaic结构的转换，坐标转换完毕后将模型导出为DXF格式。三次B样条曲线的拟合采用基于AutoCAD 2016的二次开发实现，建立3D中心线模型。已完成Matlab与二次开发相关代码编写。

(4)血管模型设计

本发明依据统计出的截面特征，建立的中心线模型，采用AutoCAD的二次开发对3D中心线模型进行批量扫掠与抽壳操作，二次开发采用的AutoLisp语言，已完成相关代码编写。

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。