三维主动脉瓣膜的体参数化构建

期刊菜单

三维主动脉瓣膜的体参数化构建
Constructing Volumetric Parametrization Model of 3D Aortic Valve

DOI: 10.12677/mos.2025.143238, PDF, HTML, XML,
作者: 刘亮：上海理工大学机械工程学院，上海
关键词: 主动脉瓣膜；等几何分析；体参数化模型；非均质材料；Aortic Valve； Isometric Analysis； Volumetric Parameterization Model； Heterogeneous Material

摘要: 主动脉瓣膜是由复杂曲面构成的，其几何模型需要高阶的自由曲面表示，其物理场需要网格离散表示，二者存在表示格式不兼容的问题，针对该问题，本文提出一种可用于等几何分析的三维主动脉瓣膜体参数化构建方法。首先，从CT图像中分割出主动脉瓣膜区域并转化为点云数据，将点云数据由十二条切割线分割为六个表面区域；其次，使用样条曲面构建每个表面区域的几何表达，由表面样条曲面经过体布尔和运算生成体参数化模型；最后，由点云数据的Hu (Hounsfield Unit)值信息构建多材料分布的参数化模型。实例表明该方法能用更少自由度表达复杂曲面并储存多材料分布信息，基于体参数化模型实现主动脉瓣膜几何表达与材料分布相耦合，完成计算机辅助设计CAD与计算机辅助工程CAE集成。

Abstract: The aortic valve is composed of complex surfaces. Its geometric model requires high-order free-form surface representation, and its physical field requires mesh discrete representation. There is an incompatible representation format between the two. To address this problem, this paper proposes a three-dimensional aortic valve body parametric construction method that can be used for isogeometric analysis. First, the aortic valve area is segmented from the CT image and converted into point cloud data. The point cloud data is segmented into six surface areas by twelve cutting lines; secondly, the spline surface is used to construct the geometric expression of each surface area, and the surface spline surface is used to generate a volume parametric model through volume Boolean and operation; finally, the Hu (Hounsfield Unit) value information of the point cloud data is used to construct a multi-material distribution parametric model. Examples show that this method can express complex surfaces with fewer degrees of freedom and store multi-material distribution information. Based on the volume parametric model, the aortic valve geometric expression and material distribution are coupled, and the integration of computer-aided design CAD and computer-aided engineering CAE is completed.

文章引用：刘亮. 三维主动脉瓣膜的体参数化构建[J]. 建模与仿真, 2025, 14(3): 460-470. https://doi.org/10.12677/mos.2025.143238

1. 引言

主动脉瓣疾病不仅是老年人最常见的心脏瓣膜病，也是心脏病的第三大病因，主动脉疾病的手术治疗，不仅需要主治医生的丰富经验，也需要患者特异性的瓣膜力学性能仿真结果提供支撑。主动脉瓣膜的力学性能仿真结果的准确性在很大程度上取决于所构建模型的精确度和所选用仿真方法的可靠性。鉴于心脏瓣膜形态的复杂性，传统建模方式在构建精确模型时面临诸多挑战，尤其是在二维图像上进行疾病评估时。然而，CT影像技术可提供丰富的解剖结构、形状和尺寸数据[1]，基于这些数据，CT影像重建技术能够精确地重建主动脉瓣膜的几何形状。

在心脏瓣膜三维重构领域，我国学者已取得了显著进展。郭燕丽[2]等人于2003年首次利用Medview三维重构软件和contours + marching cubes算法，通过数据分割、光滑处理等步骤，成功实现了我国第一例人体心脏的可视化三维重构。随后，张潇[3]在2016年通过提取医学影像的扫描数据，结合阈值提取、修改编辑、区域增长及二次网格划分等方法，完成了心脏瓣膜和动脉壁的重建，并利用逆向软件和有限元分析软件对模型进行了优化和仿真分析。而陈彦飞[4]等人在2020年则通过CT影像资料详细介绍了从几何重建、逆向工程、网格划分到材料定性等建立有限元模型的详细处理步骤。利用人体三维重构应用软件Mimics [5] [6]和逆向工程软件Geomagic [7] [8]处理得到瓣膜模型，并且用于ANSYS Workbench仿真分析的方法已相对成熟。然而，尽管这些方法在构建有限元模型方面取得了显著成果，但尚未有成熟的方法能够将重建模型转化为可供等几何分析的体参数化模型。

近年来，等几何分析[9]作为一种新兴的数值分析方法，实现了CAD和CAE的一体化，有效避免了建模仿真过程中模型表达不一致的问题[10]。能够完整地保留瓣膜的形状特征，而样条技术的应用则进一步提高了患者特异性瓣膜几何形状的构建精度。

因此，本文以经CT图像构建主动脉瓣膜真实几何模型和材料信息为目标，致力于研究基于IGA的心脏瓣膜体参数化模型。总结了现有CT重建相关处理步骤及算法，并提出了一套从CT重建到兼顾几何形状和材料属性的体参数化模型完整构建流程，为实现更精确的心脏瓣膜力学性能仿真提供有力支持。

2. 点云数据重建

为了构建主动脉瓣膜几何与材料模型，本研究依据患者的临床CT图像数据实现了点云数据重建，具体步骤如下所述：

(1) 首先以医学数字成像和通信(DICOM)格式导入至图像处理软件Mimics 19.0中。鉴于心脏内部不同组织(如血液、原生瓣膜及钙化瓣膜)之间存在显著的密度差异，这些差异在CT图像上体现为各不相同的Hu值分布，从而使得心脏瓣膜的结构与形态在不同切面图像上清晰可见。图1展示了四个不同视角的图像窗口，分别为正视图(左上角)、俯视图(右上角)、侧视图(左下角)以及三维视图(右下角)，通过调整这些视图的角度以确保其与主动脉截面垂直对齐。

Figure 1. Image data import

图1. 影像数据导入

(2) 由于心脏内部不同组织器官间的密度值存在显著差异，本研究利用Mimics软件中的阈值分割功能，基于瓣膜处Hu值的差异精确提取瓣膜结构。通过不断调整阈值的最高及最低值，当像素的Hu值处于所设定的阈值范围内时，瓣膜结构将被完全选中。然而，在初步分割结果中，除了瓣膜结构外，还可能包含血液和其他非目标组织。因此，需通过手工交互删除被误选中的多余结构。经过多次重复操作后，可成功分割出主动脉瓣膜区域，如图2所示，其中黄色高亮部分即为分割结果，分别展示了瓣膜冠状面(a)、矢状面(b)和水平面(c)的视图。

Figure 2. Threshold segmentation

图2. 阈值分割

(3) 随后，采用Mimics软件中的“从蒙版计算3D”功能，将步骤(2)中分割得到的主动脉瓣膜转化为三维点云数据。然而，初步生成的三维模型表面可能会存在一些明显的钉状凸起和孔洞等缺陷，如图2所示。

3. 体参数化模型构建

首先，选中瓣膜表面点云中的数据点构造十二条切割线，并通过求解最小二乘问题得到单变量样条曲线表示的切割线。其次，利用切割线将瓣膜划分为六个表面，在每个区域上构建双变量样条曲面作为基面，分别逼近表面上的点云数据点，得到瓣膜的真实曲面。最后，基于六个表面使用体布尔和操作构建三变量样条体，得到体参数化模型。

3.1. 构造切割线

给定数据点 ${Q_{i} | i = 0, 1, \dots, n}$ ，使用 $p$ 次样条曲线插值数据点，该曲线为所有数据点的误差最小曲线

$C u r (u) = \sum_{i = 0}^{m} N_{i, p} (u) P_{i}$ ，

采用弦长参数化计算数据点 $Q_{i}$ 所对应的参数值 ${\bar{u}}_{i}$

令

$D = \sum_{i = 0}^{n} | Q_{i} - Q_{i - 1} |$ (1)

则数据点的参数值为：

${\begin{cases} u_{0} = 0 \\ u_{n} = 1 \\ u_{i} = u_{i - 1} + \frac{| Q_{i} - Q_{i - 1} |}{D}, i = 1, 2, \dots, n - 1 \end{cases}$ (2)

得到曲线的节点矢量 $U = {u_{0}, u_{1}, \dots, u_{p}, \dots, u_{n}}$ ，计算出控制点满足

${\begin{cases} Q_{(0)} = C u r (0) \\ Q_{(1)} = C u r (1) \end{cases}$ (3)

其余的控制点通过下式计算

$P = {(N^{T} N)}^{- 1} N^{T} R$ (4)

其中

$N = [\begin{matrix} N_{0, p} (u_{0}) & N_{1, p} (u_{1}) & \dots & N_{n, p} (u_{n}) \\ N_{0, p} (u_{0}) & N_{1, p} (u_{1}) & \dots & N_{n, p} (u_{n}) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ N_{0, p} (u_{0}) & N_{1, p} (u_{1}) & \dots & N_{n, p} (u_{n}) \end{matrix}], R = [\begin{matrix} Q_{1}^{x} - G^{x} & Q_{1}^{y} - G^{y} & Q_{1}^{z} - G^{z} \\ ⋮ & ⋮ & ⋮ \\ Q_{N - n}^{x} - G^{x} & Q_{N - n}^{y} - G^{y} & Q_{N - n}^{z} - G^{z} \end{matrix}]$ (5)

$G_{k} = N_{0, p} (u_{k}) Q_{0} + N_{m, p} (u_{k}) Q_{m}$ (6)

求出控制点即可得到切割线 $C u r (u)$ 。

3.2. 构造表面

一个表面上四条样条曲线 $P (u, 0)$ ， $P (u, 1)$ ， $P (0, v)$ ， $P (1, v)$ ， $u, v \in [0, 1]$ ，如图3，包围成封闭的曲边四边形，延用 $u, v$ 方向上的样条曲线节点矢量基面节点矢量 ${U}$ ， ${V}$ Coons曲面表达式为：

$\begin{array}{l} S (u, v) = (1 - v) P (u, 0) + v P (u, 1) + (1 - u) P (0, v) + u P (1, v) \\ - [(1 - v) (1 - u) P (0, 0) + (1 - v) u P (1, 0) + v (1 - u) P (0, 1) + u v P (1, 1)] \end{array}$ (7)

使用上述步骤构造的Coons曲面为投影基面，使用拟合数据点的方法，得到瓣膜六个边界面的真实曲面。

Figure 3. Quadrilateral bounding wireframe

图3. 四边形边界线框

计算出数据点在基面上的投影参数值，先采用搜索方法，将搜索范围按照参数与均匀分成两份，迭代逐步缩小搜索范围，直到数据点被包含于搜索范围内。如果采用等分节点将基面等分为 $u$ 向 $m$ 等份、 $v$ 向 $n$ 等分，则需进行 $(m + 1) (n + 1)$ 次数据点投影参数确认，若要求投影参数的精度达到0.1时 $m = 10$ ，即计算次数为121次，计算效率过低。采用折半搜索法，将参数域均匀分成两等份，一次计算后，以投影参数 $(u, v)$ 为中心，以 $[u - 1 / 2^{t}, u + 1 / 2^{t}] [v - 1 / 2^{t}, v + 1 / 2^{t}]$ 为搜索区间再将参数域分割成两份，最终计算量仅有 $6 t$ 次。确定投影的大致区域后，接着利用牛顿迭代方法精确的计算出数据点在基面投影点，当迭代后的距离与前一次相比差异小于定值，则将当前参数值视为数据点参数，并停止迭代。

定义矢量值函数：

$r (u, v) = S (u, v) - P$ (8)

相对应的标量方程：

${\begin{cases} f (u, v) = r (u, v) \cdot S_{u} (u, v) = 0 \\ g (u, v) = r (u, v) \cdot S_{v} (u, v) = 0 \end{cases}$ (9)

解标量方程，设

$δ_{i} = [\begin{matrix} Δ u \\ Δ v \end{matrix}] = [\begin{matrix} u_{i + 1} - u_{i} \\ v_{i + 1} - v_{i} \end{matrix}]$ (10)

$J_{i} = [\begin{matrix} f_{u} & f_{v} \\ g_{u} & g_{v} \end{matrix}] = [\begin{matrix} {| S_{u} |}^{2} + r S_{u u} & {| S_{u} |}^{2} + r S_{u u} f_{v} \\ S_{u} S_{v} + r S_{u u} & {| S_{v} |}^{2} + r S_{v v} \end{matrix}]$ (11)

$k_{i} = - [\begin{matrix} f (u_{i}, v_{i}) \\ g (u_{i}, v_{i}) \end{matrix}]$ (12)

在第 $i$ 步迭代时，求解线性方程：

$J_{i} δ_{i} = k_{i}$ (13)

更新参数：

${\begin{cases} u_{i + 1} = u_{i} + Δ u \\ v_{i + 1} = v_{i} + Δ v \end{cases}$ (14)

收敛准则：

${\begin{cases} | S (u_{i}, v_{i}) - p | \leq ε_{1} 判断点是否重合 \\ \frac{S^{'} \cdot (S - p)}{| S^{'} | | S - p |} \leq ε_{2} 判断余弦是否归零 \end{cases}$ (15)

使用上述收敛准则确保参数在定义域内，若超过定义域，则取近端点值。

Figure 4. Spline surface fitting of point cloud data

图4. 点云数据样条曲面拟合

数据点参数化之后，需根据最小二乘法原理，计算真实曲面如图4。已知数据点 ${Q_{l}}_{i = 1}^{N}$ 在基面的投影为 $S (u, v)$ ，寻求一张样条曲面 $S u r (u, v) = \sum_{i = 0}^{m} \sum_{j = 0}^{m} N_{i, p} (u) N_{j, q} (v) P_{i, j}$ ，使其满足：

${\begin{cases} Q_{(0, 0)} = S u r (0, 0) \\ Q_{(0, 1)} = S u r (0, 1) \\ Q_{(1, 0)} = S u r (1, 0) \\ Q_{(1, 1)} = S u r (1, 1) \end{cases}$ (16)

其余数据点 $Q_{(u_{i}, v_{i})}$ 在最小二乘原理下被逼近。即

$\min \sum_{i = 0}^{N - k} | Q_{(u_{i}, v_{i})} - S u r (u_{i}, v_{i}) |$ (17)

其中 $Q_{(u_{i}, v_{i})}$ 为数据点中第 $i$ 个点， $(u_{i}, v_{i})$ 则为其对应的投影参数。 $k$ 为曲面的角点个数。

根据B样条最小二乘法拟合数据点，构造控制点反求方程：

$P = {(N^{T} N)}^{- 1} N^{T} R$ (18)

其中：

$\begin{array}{l} N = [\begin{matrix} N_{1, p} (u_{1}) N_{0, q} (v_{1}) & \dots & N_{m - 1, p} (u_{1}) N_{0, q} (v_{1}) & \dots & N_{1, p} (u_{1}) N_{m, q} (v_{1}) & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ N_{1, p} (u_{N - n}) N_{0, q} (v_{N - n}) & \dots & N_{m - 1, p} (u_{N - n}) N_{0, q} (v_{N - n}) & \dots & N_{1, p} (u_{N - n}) N_{m, q} (v_{N - n}) & \dots \end{matrix} \\ \begin{matrix} N_{m - 1, p} (u_{1}) N_{n, q} (v_{1}) \\ ⋮ \\ N_{m - 1, p} (u_{N - n}) N_{n, q} (v_{N - n}) \end{matrix}] \end{array}$ (19)

$P = | \begin{matrix} p_{1, 0}^{x} & p_{1, 0}^{y} & p_{1, 0}^{z} \\ ⋮ & ⋮ & ⋮ \\ p_{m - 1, 0}^{x} & p_{m - 1, 0}^{y} & p_{m - 1, 0}^{z} \\ p_{0, 1}^{x} & p_{0, 1}^{y} & p_{0, 1}^{z} \\ ⋮ & ⋮ & ⋮ \\ p_{m, 1}^{x} & p_{m, 1}^{y} & p_{m, 1}^{z} \\ ⋮ & ⋮ & ⋮ \\ p_{0, n - 1}^{x} & p_{0, n - 1}^{y} & p_{0, n - 1}^{z} \\ ⋮ & ⋮ & ⋮ \\ p_{m, n - 1}^{x} & p_{m, n - 1}^{y} & p_{m, n - 1}^{z} \\ p_{1, n}^{x} & p_{1, n}^{y} & p_{1, n}^{z} \\ ⋮ & ⋮ & ⋮ \\ p_{m, n}^{x} & p_{m, n}^{y} & p_{m, n}^{z} \end{matrix} | R = [\begin{matrix} Q_{1}^{x} - G^{x} & Q_{1}^{y} - G^{y} & Q_{1}^{z} - G^{z} \\ ⋮ & ⋮ & ⋮ \\ Q_{N - n}^{x} - G^{x} & Q_{N - n}^{y} - G^{y} & Q_{N - n}^{z} - G^{z} \end{matrix}]$ (20)

其中 $G = N_{0, p} (u) N_{0, q} (v) P_{0, 0} + N_{m, p} (u) N_{0, q} (v) P_{m, 0} + N_{0, p} (u) N_{n, q} (v) P_{0, n} + N_{m, p} (u) N_{n, q} (v) P_{m, n}$ ，有 $(m + 1) (n + 1) - 4$ 个未知数 $P_{i, j}$ ， $(m + 1) (n + 1) - 4$ 个方程，可求出拟合样条曲面的控制点，控制点云数据的曲面拟合精度为数据点到面片之间的平均距离。使用上一次的拟合曲面为基面对数据点重新投影参数，再重复进行最小二乘拟合。

$e = \sqrt{\sum_{i = 0}^{n} {| Q_{i} - S u r (u_{i}, v_{i}) |}^{2} / n}$ (21)

其中， $Q_{i}$ 为数据点， $S u r (u_{i}, v_{i})$ 为重构曲面上数据点的投影点。当迭代误差差值的小于指定阈值时，将拟合曲面作为真实曲面，迭代结束。反之，若无法满足精度需求，则通过细分方法控制点个数增加，再进行最小二乘拟合迭代，直至达到精度需求。

3.3. 构造体

当六个真实表面均构造完成后，使用这六个面进行样条体控制点插值。设六个面 $S u r_{i}, i = 1, 2, \dots, 6$ ，的参数方程为：

$\begin{array}{l} U 方向： S u r_{3} = C (0, v, w), S u r_{4} = C (1, v, w) \\ V 方向： S u r_{1} = C (u, 0, w), S u r_{2} = C (u, 1, w) \\ W 方向： S u r_{5} = C (u, v, 0), S u r_{6} = C (u, v, 1) \end{array}$

六个面所对应的12条边界曲线，如图5，曲线方程分别为：

$\begin{array}{l} U 方向： L_{15} = C (u, 0, 0), L_{16} = C (u, 0, 1), L_{25} = C (u, 1, 0), L_{26} = C (u, 1, 1) \\ V 方向： L_{35} = C (0, v, 0), L_{36} = C (0, v, 1), L_{45} = C (1, v, 0), L_{46} = C (1, v, 1) \\ W 方向： L_{13} = C (0, 0, w), L_{14} = C (1, 0, w), L_{23} = C (0, 1, w), L_{24} = C (1, 1, w) \end{array}$

六个面所对应的八个顶点：

$P_{1} = C (0, 0, 0)$ ， $P_{2} = C (1, 0, 0)$ ， $P_{3} = C (1, 0, 1)$ ， $P_{4} = C (0, 0, 1)$

$P_{5} = C (0, 1, 0)$ ， $P_{6} = C (1, 1, 0)$ ， $P_{7} = C (0, 1, 0)$ ， $P_{8} = C (0, 1, 1)$

使用 $U$ 、 $V$ 、 $W$ 方向上的两曲面进行线性插值：

$C_{1} (u, v, w) = (1 - u) S u r_{3} + u S u r_{4}$

$C_{2} (u, v, w) = (1 - v) S u r_{1} + v S u r_{2}$

$C_{3} (u, v, w) = (1 - w) S u r_{5} + w S u r_{6}$

同理，由12条边界曲线线性插值得到的体模型 $U$ 、 $V$ 、 $W$ 方向上为：

$\begin{array}{l} C_{4} (u, v, w) = (1 - v) (1 - w) C (u, 0, 0) + v (1 - w) C (u, 1, 0) \\ + v w C (u, 1, 1) + (1 - v) w C (u, 1, 1) \end{array}$ (22)

$\begin{array}{l} C_{5} (u, v, w) = (1 - u) w C (0, v, 1) + (1 - u) (1 - w) C (0, v, 0) \\ + u w C (1, v, 1) + u (1 - w) C (1, v, 0) \end{array}$ (23)

$\begin{array}{l} C_{6} (u, v, w) = (1 - u) (1 - v) C (0, 0, w) + u (1 - v) C (1, 0, w) \\ + u v C (1, 1, w) + (1 - u) v C (0, 1, w) \end{array}$ (24)

八个顶点线性插值得到的体模型分别为：

$\begin{array}{l} C_{7} (u, v, w) = (1 - u) (1 - v) (1 - w) C (0, 0, 0) + u (1 - v) (1 - w) C (1, 0, 0) \\ + u (1 - v) w C (1, 0, 1) + (1 - u) (1 - v) w C (0, 0, 1) \\ + (1 - u) v (1 - w) C (0, 1, 0) + u v (1 - w) C (1, 1, 0) \\ + u v w C (1, 1, 1) + (1 - u) v w C (0, 1, 1) \end{array}$ (25)

故六面体模型的体插值公式为：

$\begin{array}{l} C (u, v, w) = C_{1} (u, v, w) + C_{2} (u, v, w) + C_{3} (u, v, w) - C_{4} (u, v, w) \\ - C_{5} (u, v, w) - C_{6} (u, v, w) + C_{7} (u, v, w) \end{array}$ (26)

Figure 5. Hexahedron boundary wireframe

图5. 六面体边界线框

4. 材料模型构建

瓣膜点云中每个数据点所携带的Hu值信息表示模型的材料特性。将离散的材料参数化，并据此构建出一个于几何模型耦合的连续材料空间。即构建与几何模型耦合的材料体参数化模型表达如下式：

$M (u, v, w) = \sum_{i = 0}^{m} \sum_{j = 0}^{n} \sum_{k = 0}^{l} M_{i j k} R_{i, p} (u) R_{j, q} (v) R_{k, r} (w)$ (27)

式中的 $R_{i, p} (u) R_{j, q} (v) R_{k, r} (w)$ 与几何模型相同，进一步的寻找控制点 $P_{i j k}^{}$ 与之欧式距离最小数据点 $Q_{l}$ 对应的Hu值 $M_{l}$ ，并将其作为 $M_{i j k}^{}$ 的初值 $M_{i j k}^{(0)}$ 。假设第 $m$ 次迭代所得的Hu值为 $M_{}^{(m)} (u_{l}, v_{l}, w_{l})$ ，则每一步迭代从 $M_{}^{(m)} (u_{l}, v_{l}, w_{l})$ 至 $M_{}^{(m + 1)} (u_{l}, v_{l}, w_{l})$ 包含以下步骤：

1) 计算数据点差分向量

$δ_{l}^{(m)} = M_{l} - M_{}^{(m)} (u_{l}, v_{l}, w_{l})$ (28)

其中参数 $(u_{l}, v_{l}, w_{l})$ 为数据点 $Q_{l}$ 至材料体参数化模型中的逆映射。

2) 每一个数据点的Hu值 $M_{l}^{(m)}$ 为相关控制点处的Hu值 $M_{i j k}^{(m)}$ 的加权组合，其组合系数为 $R_{i, p} (u_{l}) R_{j, q} (v_{l}) R_{k, r} (w_{l}) \neq 0$ ，反之， $M_{i j k}^{(m)}$ 也是相关数据点处 $M_{l}^{(m)}$ 的加权组合，因此可以得到 $M_{i j k}^{(m)}$ 的差分向量：

$Δ_{i j k}^{(m)} = \frac{\sum_{l \in I_{i j k}} R_{i, p} (u_{l}) R_{j, q} (v_{l}) R_{k, r} (w_{l}) δ_{l}^{(m)}}{\sum_{l \in I_{i j k}} R_{i, p} (u_{l}) R_{j, q} (v_{l}) R_{k, r} (w_{l})}$ (29)

其中 $I_{i j k}$ 是 $R_{i, p} (u_{l}) R_{j, q} (v_{l}) R_{k, r} (w_{l}) \neq 0$ 的集合。迭代结果的极限形式为材料模型对数据点Hu值的最小二乘拟合，当迭代矩阵的谱半径小于1时，迭代结果收敛。具体步骤可见算法1。

算法1：非均质材料模型构建

输入：带有Hu值 ${M_{l}}$ 数据点 ${P_{l}}_{l = 0, 1, \dots, m}$ ，误差 $t o l = 0.01$

输出：体参数化几何模型 $V (u, v, w)$ ，材料模型 $M (u, v, w)$

1. 初始化带有数据点集 ${{P^{'}}_{l}}_{i = 0, 1, \dots, n}$ 的曲面集 ${S_{i}^{(0)}}_{i = 0, 1, \dots, 5}$

2. 由曲面数据点及其投影点计算数据点坐标 ${ε_{l}^{(0)} = P_{l}^{(0)} - V (u_{l}, v_{l}, w_{l})}$ 的差分向量集

3. while $| {\sum_{l} ‖ ε_{l}^{(m + 1)} ‖}^{2} / {\sum_{l} ‖ ε_{l}^{(m)} ‖}^{2} | > t o l$ do

4. 计算差向量 $ε_{l}^{(m)}$ ，同式(3.31)

5. 计算 $Ε_{i j}^{(m)}$ 与 $ε_{i j}^{(m)}$ 的差向量，同式(3.32)

6. 更新每个曲面 $S_{i}$ 上的控制点 $P_{i j}^{^{(m + 1)}} = P_{i j}^{^{(m)}} + Ε_{i j}^{^{(m)}}$

7. end while

8. 使用体布尔和操作通过曲面 ${S_{i}^{(m)}}$ 构建体参数化模型 $V (ξ, η, ζ)$

9. 将控制点最近数据点的Hu值设置为控制点Hu的初始化集合 ${M_{i j k}^{^{(0)}}}$ ，同时设置数据点处的Hu值及其投影点的Hu值初始化差分向量集合 ${δ_{l}^{(0)} = M_{l} - M_{}^{(0)} (ξ_{l}, η_{l}, ζ_{l})}$

10. while $| {\sum_{l} ‖ δ_{l}^{(m + 1)} ‖}^{2} / {\sum_{l} ‖ δ_{l}^{(m)} ‖}^{2} | > t o l$ do

11. 通过式(3.31)计算差向量 $δ_{l}^{(m)}$

12. 通过式(3.32)计算 $Δ_{i j}^{(m)}$ 与 $δ_{l}^{(m)}$ 的差向量

13. 更新控制点上的Hu值 $M_{i j k}^{^{(m + 1)}} = M_{i j k}^{^{(m)}} + Δ_{i j k}^{(m)}$

14. end while

15. 使用式(3.30)计算材料模型 $M (u, v, w)$

16. return体参数化几何模型 $V (u, v, w)$ ，材料模型 $M (u, v, w)$

5. 算法实例

预处理步骤中，将CT图像中瓣膜的区域分割后，根据每一个像素点的坐标建立起点云数据。构建体参数化模型时，首先构建由点云数据中的切割线，十二条切割线将瓣膜表面点云数据分割为六个区域，接着在每个表面区域上通过边界构造基面，计算点云数据在基面上的投影参数，用最小二乘原理拟合投影参数得到真实曲面，然后由真实曲面构造体模型，使用体布尔和运算构建三变量样条体，即对瓣膜进行提参数化，最后通过点云的Hu值构建瓣膜参数化的材料属性，得到多材料分布的体参数化模型。

Figure 6. Volume parameterization model construction process

图6. 体参数化模型构建过程

图6(a)是切割线生成后的点云数据，切割线将瓣膜的点云数据分割为拓扑六面体，使其满足构建三变量样条体的要求。图6(b)是经过拟合后生成的六个区域的真实曲面可以看出其精准的表示了复杂曲面。图6(c)是通过真实曲面构建的体参数化模型，根据等参线的分布可以看出其具有不错的模型质量，可用于等几何分析。图6(d)是真实材料属性分布图，根据颜色条的分布可以看出瓣膜的集中钙化区域(红色)。

Figure 7. Finite element model construction process

图7. 有限元模型构建过程

图7(a)是通过点云模型构建的三角网格构建模型，可以看出三角网格十分密集程度。图7(b)是通过三角网格构建的四面体网格有限元模型，其节点数量和计算单元数量的对比在表1。

Table 1. Comparison of volumetric parametrization model and finite element model

表1. 体参数化模型与有限元模型的对比

模型类型	节点(控制点)数量	计算单元数量
体参数模型	500	192
实体网格	8176	126,499

6. 结论

本文提出一种主动脉瓣膜体参数化模型的构建方法。通过CT图像提取出点云数据、构建分割线、真实曲面以及三变量样条体，同时构建与几何模型耦合的多材料分布参数化信息。此方法可完成表面复杂的心脏瓣膜曲面模型体参数化构建及多材料分布表达，实现CAD建模和CAE分析的一体化集成。然而，模型分割中的手工交互过程在一定程度上影响了构建效率，未来需要智能方法加以改善。

参考文献

[1]	董双国. 浅谈主动脉夹层的CT影像诊断[J]. 中国医学工程, 2012, 20(2): 173-174.
[2]	郭燕丽, 张绍祥, 刘正津, 等. 首例中国可视化人体心脏三维重建及临床意义[J]. 第三军医大学学报, 2003, 25(7): 569-571.
[3]	张潇. 基于CT图像的生物瓣膜三维模型重建及有限元分析[D]: [硕士学位论文]. 济南: 山东大学, 2016.
[4]	陈彦飞, 鲁超, 赵勇, 等. 基于CT影像动态膝关节有限元模型的构建及仿真力学分析[J]. 中国骨伤, 2020, 33(5): 479-484.
[5]	王娇, 刘洋, 张晓玲, 等. Mimics软件在医学图像三维重建中的应用[J]. 医疗卫生装备, 2015, 36(2): 115-118.
[6]	张海峰, 史冠忠, 朱洪智, 等. 基于图像数据建立器官三维模型[J]. 中国组织工程研究, 2015, 19(40): 6547-6552.
[7]	阚晓昕, 白一帆, 马韬, 等. 基于CT影像重建模型的主动脉夹层流固耦合模拟[J]. 水动力学研究与进展: A辑, 2015, 30(6): 701-706.
[8]	Li, J., Yan, W., Wang, W., Wang, S. and Wei, L. (2022) Comparison of Balloon-Expandable Valve and Self-Expandable Valve in Transcatheter Aortic Valve Replacement: A Patient-Specific Numerical Study. Journal of Biomechanical Engineering, 144, Article ID: 104501. https://doi.org/10.1115/1.4054332
[9]	Hughes, T.J.R., Cottrell, J.A. and Bazilevs, Y. (2005) Isogeometric Analysis: CAD, Finite Elements, NURBS, Exact Geometry and Mesh Refinement. Computer Methods in Applied Mechanics and Engineering, 194, 4135-4195. https://doi.org/10.1016/j.cma.2004.10.008
[10]	陈龙, 张军洋, 李郝林, 等. 面向材料反求的非均质体参数化模型构建[J]. 机械工程学报, 2014, 50(23): 156-164.

为你推荐

友情链接