基于显微CT测量的缺损颅骨组织工程血管支架参数化设计方法与流程

本发明属于组织工程技术领域，具体涉及一种基于显微CT测量的颅骨组织工程血管支架参数化设计方法。

背景技术：

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

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

技术实现要素：

为了解决现有技术中的问题，本发明提供一种基于显微CT测量的颅骨组织工程血管支架参数化设计方法，解决现有技术中颅骨组织工程支架由于缺乏血管结构从而限制了细胞的生长及成骨量的形成的问题。

本发明的技术方案是：

一种基于显微CT测量的缺损颅骨组织工程血管支架参数化设计方法，包括以下步骤：

(1)获取天然颅骨显微CT数据并进行三维重建；

利用三维可视化软件Avizo对图像进行三维重建，获得天然颅骨样本三维模型，然后对重建的体数据进行处理，只表达板障中的空腔结构；对得到的板障空腔结构，从中区分出主静脉和小静脉结构；

(2)建立血管支架参数化模型；

定义血管网络由若干条血管段相互连接而成，每条血管段可看作是由径向截面沿轴向中心线扫掠而成的曲面，通过对主静脉中心线以及径向截面进行提取与参数分析，据此构建颅骨血管的参数化模型，包括中心线参数化模型、血管段参数化模型和参数范围；

(3)颅骨血管支架参数化设计；

基于步骤(2)参数化模型，首先应用Matlab生成血管二叉树的分叉点，再运用MFC/ACIS/HOOPS相结合的图形开发平台插值拟合血管中心线段和构建扫掠截面，最后通过UG软件进行曲面缝补、加厚和在分叉处倒圆角处理，可实现血管支架三维参数化设计。

所述步骤(2)中中心线参数化模型用二叉树模型来描述中心线分叉点的分布特征，并用这些分叉特征点拟合出血管中心线曲线，二叉树节点间的位置关系定义为：

二叉树从根节点，即第0级节点开始，向上逐级分叉，通过以下参数的控制，在特定空间中可生成血管的分叉特征点；定义空间X为血管支架的设计范围尺寸，在空间X的中心处建立坐标系O-XYZ，假定二叉树第i级节点为P_i，其父节点为P_i-1，两个子节点分别为这四个点则构成基本的分叉单元，基于每个基本分叉单元，由父节点逐级生成子节点，其中每个父节点最多有两个子节点，没有子节点的父节点为端节点，这样则可生成整个二叉树模型；

定义P_i-1P_i为主支，和为分支，过两个分支的平面定义为分支平面，XOY平面为中间分层面，分支平面与中间分层面之间若有较小的夹角，可确保分叉点分布在中间分层面附近，以提高分叉级别；为空间X边界上的点，为在空间X的最大延伸长度，较大的空间延伸长度，可确保分叉点在空间X中均匀填充；节点间距离为d，θ为主支与分支间的夹角，α为分支之间的夹角，即为向量P_i-1P_i与向量的夹角，α_i为向量与向量之间的夹角；节点的坐标与夹角之间满足如下关系；

计算每个基本分叉单元的四点坐标值，由上述公式可得到主支与分支之间的夹角θ、两个分支之间的夹角α以及节点距d的参数范围。

所述步骤(2)中中心线曲线参数化模：依据基本分叉单元及参数条件可在空间X中随机生成血管分叉点，假定所有分叉点共确定了N条血管段，每条血管段有n个特征点(P₀ P₁P₂ P₃ … P_n-2 P_n-1)，定义向量V_s为中心线的起始点P₀处的切矢，向量V_e为末端点P_n-1处的切矢，且：

V_s＝P₀P₁ (4)

V_e＝P_n-2P_n-1 (5)

所述步骤(2)中血管段曲面扫掠参数化模型：定义血管半径为r，整体血管模型的初始半径为R₀，令R₀＝r_max；对于每条血管段，在初始点处确定初始扫掠截面，沿着中心线以一定锥度扫掠，可得半径逐渐变化的血管段扫掠曲面；扫掠截面圆心为P₀，法矢为V_s，扫掠锥度为γ，初算tanγ＝(r_max-r_min)/l'，l'为中心线模型的最长长度，扫掠半径为R，每条血管段在第i个特征点处

r_i＝R-l·tanγ (6)

式中，R为某一条血管段的起始半径，每条血管段的R可依据R₀逐级计算得到；l为起始点到第i个特征点间的长度；d_i为节点段长度；据此，可构建血管段扫掠曲面。

所述步骤(2)中参数范围为：

所述步骤(3)中颅骨血管支架参数化设计的实现过程为：

(1)生成二叉树模型

血管二叉树中分支点的位置确定主要包括以下步骤：

a)生成随机点只在X的边界面附近生成点距为0.1mm的随机点，分叉点与每个随机点的连线方向都代表一种可能的分支方向，由于随机点在空间的边界面附近，可保证分支点在其分支方向上有较大的空间延伸长度，从而使在颅骨支架空间内均分布有血管；

b)筛选分支方向基于基本分叉单元，由父节点同时筛选出满足θ和α参数范围的两个分支方向，并从中选择分支平面与中间分层面夹角最小的一组作为最终的分支方向；

c)确定分支点位置在筛选到的分支方向上，根据d的范围，随机生成两个分支长度，确定两个子节点的最终位置；

d)干涉判断空间生成的中心线可能存在相交，因此需要进行干涉判断；如果两条节点段之间的最短距离小于这两条血管半径之和，则生成的血管段曲面会发生干涉；为避免干涉，首先依据R₀设定阈值ε，然后计算新生成的节点段与已生成的节点段之间的最短距离，比较该距离与ε的大小，如果小于ε，则认为该分支点为干涉点，将其删除；

e)重复(b)～(d)步骤，逐级生成两个子节点，直至在X中最大限度地生成二叉树结构；

(2)拟合血管段中心线曲线

分叉点以二叉树形式生成之后，需将节点分配到血管段中，在此通过遍历二叉树方法实现：首先定义搜索方向，二叉树从低级到高级搜索为正向搜索，反之为逆向；然后正向搜索端节点，作为每条血管段末端节点；最后逆向逐级搜索父节点，直到该节点与其它血管段中的节点重复，则完成一条血管段的节点分配；

依据公式(4)、(5)计算每条血管段在P₀，P_n-1处的切矢V_s，V_e，假定每段血管段上分配有n个节点，利用三次样条曲线模型将n个型值点拟合成n-1段曲线段，且各段曲线段在型值点处要满足二阶连续的条件，第j段曲线段的数学表达式为

Q_j(t)＝T·M_h·G_h,t∈(0,1) (9)

式中，T＝(t³ t² t 1)；为Hermite矩阵；

G_h＝(P_j-1 P_j R_j-1 R_j)^T,P_j-1和P_j为型值点坐标值，R_j-1、R_j为P_j-1、P_j处的切矢；据此拟合血管中心线曲线；

(3)扫掠血管曲面

曲面参数化模型，求取初始截面半径R，法矢V_s，利用造型内核ACIS的圆形平面构造函数api_make_planar_disk，在初始节点处创建扫掠截面；然后利用扫掠函数api_sweep_with_options，以锥度γ沿血管段中心线轴向扫掠，生成主血管初始实体模型；最后利用UG，将所有血管段扫掠曲面缝合为连续曲面，加厚连续曲面，并在分叉处进行倒圆角处理，实现主血管的最终设计。

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

附图说明

图1颅骨板障静脉解剖学示意图，A、B为人体板障静脉中两种不同的形态；

图2(a)是天然颅骨扫描样本，图2(b)是重建样本体模型，图2(c)是板障中的空腔结构；

图3(a)是重建后的主静脉，图3(b)是光顺后的主静脉，图3(c)是分叉点提取示意图；

图4是分叉单元示意图；

图5(a)是A段静脉截面曲线提取示意图，图5(b)是B段静脉截面曲线提取示意图；

图6(a)是中心线曲线参数化模型，图6(b)是锥度示意图，图6(c)是血管曲面扫掠结果；

图7(a)是主血管中心线拟合结果，图7(b)是主血管扫掠结果，图7(c)是主血管与小血管中心线网络，图7(d)是最终血管模型；

图8(a)是实例一血管中心线网络，图8(b)是实例一血管设计模型；

图9(a)是实例二血管中心线网络，图9(b)是实例二血管设计模型。

具体实施方式

下面结合附图对本发明进行详细说明。

本发明的基于显微CT测量的缺损颅骨组织工程血管支架参数化设计方法，包括以下步骤：

(1)获取天然颅骨显微CT数据并进行三维重建

选取人体自然颅骨额骨部位(由天津医科大学提供)，尺寸约100mm×70mm×4mm，扫描区域如图2(a)所示。通过显微CT测量(天津三英精密仪器有限公司所提供的X射线三维显微镜nanoVoxel设备)，得到图像共计720帧，其图像分辨率为49.8μm。运用利用三维可视化软件Avizo对图像进行三维重建，获得天然颅骨样本三维模型如图2(b)所示，为了对内部的板障静脉特征进行识别，先对重建的体数据进行均匀背景和滤波处理，去除部分噪声和伪影，然后用阈值分割和不透明度处理，将表征噪声和骨组织(即颅骨内板、板障和外板)的灰度区域设为透明，只表达板障中的空腔结构，结果如图2(c)所示。由图可直观看出，该空腔结构主要是骨组织间的孔洞、管径较粗的主静脉和管径较细的小静脉。

对得到的板障空腔结构，需要从中区分出主静脉和小静脉结构，这里运用Avizo软件，通过提取其“骨架”，可识别出板障静脉。骨架是指图像中目标物体的骨骼部分，可以表征物体的拓扑结构。应用Avizo的体元细化算法，将图2(c)模型划分为若干个微小段，每个微小段作为一个体元，每个体元可提取一条骨架线，该骨架线上的体素到该体元边界的最短距离相等，该距离即每个体元的细化半径，对细化半径范围统计分析，得到主静脉管径范围为0.8～1.7mm；小静脉管径范围为0.5～0.8mm。

(2)建立血管支架参数化模型

通过图2(c)可以看出板障静脉是由主静脉和小静脉连接而成的血管网络，且小静脉与主静脉有相似的拓扑结构。以主静脉为例，如果定义血管网络由若干条血管段相互连接而成，每条血管段可看作是由径向截面沿轴向中心线扫掠而成的曲面，那么通过对主静脉中心线以及径向截面进行提取与参数分析，即可据此构建颅骨血管的参数化模型。

1)中心线参数化模型

血管网络的分叉点既反映血管整体的分布规律，又可作为每条血管段的特征点，而主静脉孔径较大，连通性较好，为此，重建主静脉曲面模型做中心线及分叉点特征分析。运用逆向工程软件Geomagic进行重建，结果如图3(a)所示，对其表面进行光顺处理，得到图3(b)所示的模型，提取其分叉点，并将分叉点用直线相连，结果如图3(c)所示。分析发现，主静脉的分叉点多为二分叉或三分叉结构，分叉点之间的距离和连线段之间的夹角，反映了分叉点在空间中的位置关系，因此可用二叉树模型来描述中心线分叉点的分布特征，并用这些分叉特征点拟合出血管中心线曲线，而中心线由节点相连而成，因此需要定义二叉树节点间的位置关系。

二叉树从根节点，即第0级节点开始，向上逐级分叉，通过以下参数的控制，在特定空间中可生成血管的分叉特征点。定义空间X为血管支架的设计范围尺寸，在空间X的中心处建立坐标系O-XYZ，假定二叉树第i级节点为P_i，其父节点为P_i-1，两个子节点分别为这四个点则构成基本的分叉单元，如图4所示。基于每个基本分叉单元，由父节点逐级生成子节点，其中每个父节点最多有两个子节点，没有子节点的父节点为端节点，这样则可生成整个二叉树模型。

图中定义P_i-1P_i为主支，和为分支，过两个分支的平面定义为分支平面，XOY平面为中间分层面，分支平面与中间分层面之间若有较小的夹角，可确保分叉点分布在中间分层面附近，以提高分叉级别。为空间X边界上的点，为在空间X的最大延伸长度，较大的空间延伸长度，可确保分叉点在空间X中均匀填充。节点间距离为d，θ为主支与分支间的夹角，α为分支之间的夹角，即为向量P_i-1P_i与向量的夹角，α_i为向量与向量之间的夹角。节点的坐标与夹角之间满足如下关系。

计算每个基本分叉单元的四点坐标值，由上述公式可得到主支与分支之间的夹角θ、两个分支之间的夹角α以及节点距d的参数范围。

2)血管段参数化模型

以图3(b)中椭圆标注的A、B两段主静脉为例，分别用5个等间距平面与静脉曲面相交，得到5个截面的轮廓曲线，对截面轮廓曲线进行光顺处理，发现从下往上沿着轴向方向截面有逐渐减小的趋势，为简化参数化造型，将每条血管截面看作是由若干半径变化的圆构成，如图5(a)、(b)所示。血管曲面模型是由若干血管段的扫掠截面沿中心线曲线向轴向扫掠得到，下面分别阐述中心线曲线模型和血管段曲面扫掠模型。

中心线曲线参数化模：依据基本分叉单元及参数条件可在空间X中随机生成血管分叉点，假定所有分叉点共确定了N条血管段，每条血管段有n个特征点(P₀ P₁ P₂ P₃ … P_n-2 P_n-1)，定义向量V_s为中心线的起始点P₀处的切矢，向量V_e为末端点P_n-1处的切矢，且：

V_s＝P₀P₁ (4)

V_e＝P_n-2P_n-1 (5)

血管段曲面扫掠参数化模型：定义血管半径为r，r_min～r_max的取值范围见下文，整体血管模型的初始半径为R₀，令R₀＝r_max。对于每条血管段，在初始点处确定初始扫掠截面，沿着中心线以一定锥度扫掠，可得半径逐渐变化的血管段扫掠曲面。如图6(a)、(b)所示，扫掠截面圆心为P₀，法矢为V_s，扫掠锥度为γ，初算tanγ＝(r_max-r_min)/l'，l'为中心线模型的最长长度，扫掠半径为R，每条血管段在第i个特征点处

r_i＝R-l·tanγ (6)

式中，R为某一条血管段的起始半径，每条血管段的R可依据R₀逐级计算得到；l为起始点到第i个特征点间的长度；d_i为节点段长度。据此，可构建血管段扫掠曲面，结果如图6(c)所示。

3)参数范围

综上，若给定血管壁厚h，则主血管参数化模型(Vessel Model，VM)可用VM(X,d,θ,α,r,γ,h)表示。由于主静脉与小静脉有相似的拓扑结构，该参数化模型同样适用于小血管，则小血管模型VM'可用VM'(X,d',θ,α,r',γ,h)表示。

依据临床缺损颅骨的面积大小以及颅骨板障厚度可确定X范围，大约为10mm×10mm×3mm～80mm×80mm×3mm；θ、α以及d的参数范围可由式(1)～(3)计算获得；主血管半径r与小血管半径r'参考上述统计的半径范围；γ依据公式(6)、(7)估算；血管的壁厚h与管径之比为相对厚度，其满足公式(8)。

h/2r＝0.04～0.1 (8)

式中，r为血管半径，据此可求取h范围。

综上，可得血管设计模型的参数取值范围，如表1所示。

表1血管设计参数取值范围

(3)颅骨血管支架参数化设计

基于上述参数化模型，首先应用Matlab生成血管二叉树的分叉点，再运用MFC/ACIS/HOOPS相结合的图形开发平台插值拟合血管中心线段和构建扫掠截面，最后通过UG软件进行曲面缝补、加厚和在分叉处倒圆角处理，可实现血管支架三维参数化设计。以下同样以主血管为例，说明设计的实现过程。

1)生成二叉树模型

血管二叉树中分支点的位置确定主要包括以下步骤：

a)生成随机点只在X的边界面附近生成点距为0.1mm的随机点，分叉点与每个随机点的连线方向都代表一种可能的分支方向，由于随机点在空间的边界面附近，可保证分支点在其分支方向上有较大的空间延伸长度，从而使在颅骨支架空间内均分布有血管；

b)筛选分支方向基于基本分叉单元，由父节点同时筛选出满足θ和α参数范围的两个分支方向，并从中选择分支平面与中间分层面夹角最小的一组作为最终的分支方向；

c)确定分支点位置在筛选到的分支方向上，根据d的范围，随机生成两个分支长度，确定两个子节点的最终位置；

d)干涉判断空间生成的中心线可能存在相交，因此需要进行干涉判断。如果两条节点段之间的最短距离小于这两条血管半径之和，则生成的血管段曲面会发生干涉。为避免干涉，首先依据R₀设定阈值ε，然后计算新生成的节点段与已生成的节点段之间的最短距离，比较该距离与ε的大小，如果小于ε，则认为该分支点为干涉点，将其删除；

e)重复(2)～(4)步骤，逐级生成两个子节点，直至在X中最大限度地生成二叉树结构。

2)拟合血管段中心线曲线

分叉点以二叉树形式生成之后，需将节点分配到血管段中，在此通过遍历二叉树方法实现：首先定义搜索方向，二叉树从低级到高级搜索为正向搜索，反之为逆向；然后正向搜索端节点，作为每条血管段末端节点；最后逆向逐级搜索父节点，直到该节点与其它血管段中的节点重复，则完成一条血管段的节点分配。表2反映了节点的搜索方向及分配示意图。

表2节点分配示意

依据公式(4)、(5)计算每条血管段在P₀，P_n-1处的切矢V_s，V_e，假定每段血管段上分配有n个节点，利用三次样条曲线模型将n个型值点拟合成n-1段曲线段，且各段曲线段在型值点处要满足二阶连续的条件，第j段曲线段的数学表达式为

Q_j(t)＝T·M_h·G_h,t∈(0,1) (9)

式中，T＝(t³ t² t 1)；为Hermite矩阵；

G_h＝(P_j-1 P_j R_j-1 R_j)^T,P_j-1和P_j为型值点坐标值，R_j-1、R_j为P_j-1、P_j处的切矢。据此拟合血管中心线曲线，结果如图7(a)所示。

3)扫掠血管曲面

曲面参数化模型，求取初始截面半径R，法矢V_s，利用造型内核ACIS的圆形平面构造函数api_make_planar_disk，在初始节点处创建扫掠截面；然后利用扫掠函数api_sweep_with_options，以锥度γ沿血管段中心线轴向扫掠，生成主血管初始实体模型；最后利用UG，将所有血管段扫掠曲面缝合为连续曲面，加厚连续曲面，并在分叉处进行倒圆角处理，实现主血管的最终设计，如图7(b)为所示。

用同样的方法，在主血管模型中可进一步生成小血管结构，具体过程在此不再赘述，如图7(c)和图7(d)所示。小血管结构可增加血管密度以及单位空间的血管体积，并在形态上与天然板障静脉更相似，最终生成的包括主血管和小血管的血管网模型可与颅骨多孔支架进行复合，使支架内具有血管结构。

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

实施例1参数取值如表3所示，设计结果如图8所示。该模型共有111条血管段构成，主血管体积为461.63mm³，血管总体积为928.33mm³，该空间内血管密度为8.59％。

表3实例一血管模型设计参数

实施例2为了表明节点段长度d与主支分支夹角θ对二叉树节点分布形态的不同影响，在合理范围内缩小θ和d的范围，如表4所示，设计结果如图9所示。实例二共有136条血管段构成，主血管体积为425.10mm³，血管总体积为921.27mm³，空间内血管密度为8.53％。

表4实例二血管模型设计参数

由上述两实例可看出，两组参数所得到的主血管体积、总体积以及血管密度均相差不大，主血管在空间中分布都较为均匀，小血管也都能充分地填充主血管之间的空隙，有效地提高了空间中的血管密度。缩小θ和d的范围，会缩短单条血管段长度，增加分叉点数量，可使血管分布更为紧密。由此可见，该参数化设计方法能使血管均匀、充分地分布在设定空间中，作为颅骨血管支架是可行的。

本发明的血管支架参数化设计已编写了完整的Matlab程序，根据所设计的支架尺寸要求只需输入血管支架的8个参数(X,d,d',θ,α,r,r',γ)，即可确定主血管以及小血管所有分叉点位置，并将分叉点分配到血管段中，然后进一步计算血管段数量N和每条血管段上分叉点数量n，分叉点坐标P，起始截面半径R，起始向量V_s和末端向量V_e，并把这些参数保存在文本文档中，用来进行后续血管段中心线拟合和曲面扫掠。本发明的血管支架参数化设计已建立了MFC/ACIS/HOOPS相结合的图形开发平台，只需依次读取保存在文本文档中的相关参数，即可生成基于三次样条曲线拟合的血管段中心线，并沿中心线扫掠血管段曲面，生成整个血管网络的曲面模型。利用UG在每个分叉处进行倒圆角，将曲面缝合成连续曲面，并基于壁厚h，将曲面加厚为实体模型。

以上所述仅为本发明的较佳实施例，本发明的保护范围不限于上述的实施案例，凡依本发明申请专利范围所做的均等变化和修饰，皆应属本发明的涵盖范围，本申请所要求的保护范围如本申请权利要求书所示。