0 引言

快速成型技术也称为3D打印技术，是20世纪80年代末兴起的一门新技术，近年来发展十分迅猛，已被视为一项颠覆传统生产方式的革命性技术。随着3D打印技术和医疗技术的深度结合与发展，3D打印作为一种快速原型技术，其最大优势在于对复杂结构的一体化制造^[1-2]，并且可以实现针对特定患者、特定需求的各种器官的个性化生产。由于骨组织结构相对比较简单，近年来，3D打印在骨组织修复方面的应用较多，综述性文章也较多。王镓垠等^[3]对人体组织3D打印的最新进展作了综述，主要介绍了3D打印在骨组织、血管和人工肝脏的最新研究成果。Fedorovich等^[4]以及Bose等^[5-6]论述了在3D打印骨组织工程支架的最新进展以及当前的挑战和未来的方向。Jariwala等^[7]也作了类似的综述。这些综述重点讨论了人工骨材料的材料应用以及人工骨3D打印工程的前景和挑战，并没有涉及到3D打印骨组织的实现过程。早在2000年，颜永年等^[8]就提到了人工骨的快速成型制造技术的流程：由数字成像技术获得的层片数字图像, 经处理后获得层片轮廓信息，建立三维几何模型，再进行快速成型。江静等^[9]阐述了逆向工程及快速成型技术在面容多发性骨折缺损骨和下颌骨缺损三维模型重建等医学上的应用。吴纪楠等^[10]通过快速成型技术实现了个体化下颌人工骨的临床应用，达到有效的临床效果。Kouhi等^[11]采用熔融层积成型技术打印患者下颌骨模型，也取得了不错的临床反馈。

文献[9-11]虽然都提到了结合医学图像三维重建与快速成型技术结合实现医疗快速成型的方法，并且临床实验的反馈表明人工骨可以有效改善手术规划提高准确性，但是并没有给出直接从二维医学图像序列生成模型轮廓的具体细节，而是主要依赖于相关的专业软件来实现人工骨组织的3D打印。为了探索从医学图像序列打印人体原型的过程中的关键技术的细节，本文提出了一种可以直接从医学图像序列生成骨组织模型轮廓线的简化三角片模型分层过程的方法。首先, 通过分析医学图像面绘制模型与三角网格文件 (STereoLithography, STL) 模型的本质关系，根据医学图像重建移动立方块 (Marching Cube, MC) 算法^[12]中三角片建立的顺序对三角片集合分组，分组后减少了切平面对模型遍历检索的三角片数量，可以提高切片分层速度; 然后, 采用对边追踪的方法计算每一个切平面与其对应的三角片数组的交点轮廓数据。本文采用该方法对三组骨组织图像序列进行分层处理，并分别与其对应的医学图像数据生成的STL文件的分层^[13]时间进行比较，进行了相关实验; 最后, 通过对这些轮廓数据进行填充后在3D打印机上分别打印出相应的骨骼模型。

1 MC面绘制顺序重建三角片模型

面绘制是基于图形学的绘制方式，对形体的表示、操作和显示都是基于图形学的基本元素:点、线、面以及法向量来完成的。移动立方块 (MC) 算法是一种经典的面绘制重建算法，已经得到了广泛的应用。MC方法本质上是从三维数据场中提取一个等值面^[14]，该算法的核心就是根据设定的阈值，从给定的断层图像序列中提取出等值面的三角片网格模型。

二维断层图像顺序读入后，每相邻两幅切片图像对应的四个点构成一个立方体，如图 1所示。其中左图表示二维图像序列组成的三维离散数据场，右图表示离散数据场中8个相邻的数据点构成的最小立方体，对每个立方体的8个顶点和12条边分别进行编号。

通过对每个立方体的8个顶点的灰度值与设定的等值面阈值进行比较，可以判断出等值面是否和该立方体相交，即如果立方体一条边的两个顶点的灰度值分别大于、小于给定的等值面阈值，则在这条边上有且仅有一点与等值面相交，该点为等值点。假设一个顶点的灰度值大于阈值，则将它标记；如果小于阈值，则不标记。因此，一个立方体的8个顶点分别有标记和不标记两种状态，从而建立一个8位二进制索引值，如图 2所示。

确定立方体内三角片的组成方式后，计算三角片的顶点和法向量。设p₁和p₂是立方体上与等值面相交边的两顶点的坐标值，其灰度值分别为I₁和I₂, 法向量分别为n₁和n₂，等值面阈值为T。可以采用线性插值法计算交点p及其法向量n，其计算公式分别为式 (1) 和式 (2) :

$p={{p}_{1}}+(T-{{I}_{1}})({{p}_{2}}-{{p}_{1}})/({{I}_{2}}-{{I}_{1}})$

(1)

$\mathit{\boldsymbol{n}}={{\mathit{\boldsymbol{n}}}_{1}}+(T-{{I}_{1}})({{\mathit{\boldsymbol{n}}}_{2}}-{{\mathit{\boldsymbol{n}}}_{1}})/({{I}_{2}}-{{I}_{1}})$

(2)

其中，在整个断层图像序列构成的矩阵I中，立方体某顶点q(i, j, k) 处的单位法向量n_q可由该点梯度g=(g_i, g_j, g_k) 计算得到。顶点q处的梯度可以通过该点处灰度差分计算得到，如式 (3) :

$\left\{ \begin{align} & {{g}_{i}}=[I(i+1,j,k)-I(i-1,j,k)]/2 \\ & {{g}_{j}}=[I(i,j+1,k)-I(i,j-1,k)]/2 \\ & {{g}_{k}}=[I(i,j,k+1)-I(i,j,k-1)]/2 \\ \end{align} \right.$

(3)

通过上述MC算法算理可以将二维图像序列中构成等值面的三角片模型的顶点坐标及其法向量求解出来，以用来求解重建后模型的分层轨迹。

2 医学图像重建模型分层轨迹计算

医学图像重建模型分层算法的基本思想是根据图像序列数据中最小立方体的顺序建立三角片的分组矩阵，从而有效地减少切平面对模型遍历检索的三角片数量，提高切片分层速度。最后通过对边追踪^[15-16]的方法对分组后的三角片矩阵按照顺序求解切平面对边相交的轮廓线。

2.1 三角片分组及三角片表建立

在医学图像序列面绘制重建过程中，由于是对图像数据顺序处理生成的三角片模型，每一层的图像数据对应一个平面高度，假设切平面从下往上对模型分层，则与三角片模型相交的切片平面出现的次序也不同, 位置低的三角片先被切到, 位置高的后被切到。因此可以直接对三角片生成的顺序分组，并不占用额外的排序时间，即每一层图像数据生成的三角片集合对应一个数组。三角片数据整体分组的过程如图 3所示。

在图 3左侧中，从下往上是按照图像数据顺序生成的面绘制模型，依次在Z方向上增加三角面片，每一层图像数据对应一个分组序号t，每一个序号对应的三角片集合为{t_nm1, t_nm2, …, t_nmn}，如图 3右侧所示。整体分组后，假设每组三角片集合在Z方向的最小值为z_min，最大值为z_max，当某切平面的高度z_i满足z_min＜z_i＜z_max时，则该切平面仅与这一组三角片集合有交点，不需要检索其他组内三角片，并且当分层间距设定的比较大时，有些三角片在分层切片时并不参与计算，可以直接省去，大大减少了三角片的检索数量。

2.2 切平面与对应三角片数组求交

建立好三角片分组数组之后，开始计算每一个切平面与其对应的三角片相交的轮廓线。采用对边追踪的方法无需预先建立STL模型所有三角片的拓扑面关系。求解过程中, 对于一个切片平面z_i, 首先找到第1个与该切片平面相交的三角片t_ij，判断出该三角片三个顶点中在Z方向上的值最小的顶点和值最大的顶点，如果该三角片Z方向上最大值和最小值相同并且等于切平面z_i的高度zValue，为避免过多冗余点产生，则忽略该三角片，继续其他三角片的计算；如果不相等，则可以通过式 (4) 计算切平面与三角片的第一个交点坐标p₁(x_p, y_p, z_p)：

$\left\{ \begin{align} & {{x}_{p}}=\frac{{{z}_{p}}-{{z}_{\min }}}{{{z}_{\max }}-{{z}_{\min }}}\left( {{x}_{\max }}-{{x}_{\min }} \right)+{{x}_{\min }} \\ & \\ & {{y}_{p}}=\frac{{{z}_{p}}-{{z}_{\min }}}{{{z}_{\max }}-{{z}_{\min }}}\left( {{y}_{\max }}-{{y}_{\min }} \right)+{{y}_{\min }} \\ & \\ & {{z}_{p}}=zValue \\ \end{align} \right.$

(4)

其中：x_min、y_min、z_min和x_max、y_max、z_max分别代表最小顶点p_Min的坐标和z值最大顶点p_Max的坐标；起始时，p_Min代表起始三角片z值最小的顶点，p_Max代表起始三角片z值最大的顶点；随后最小顶点p_Min和最大顶点p_Max则根据对边追踪的结果不断改变，以便计算切平面与相邻三角片的交点坐标，追踪过程如下：

1) 检索分组三角片矩阵中的三角片顶点坐标。

2) 该三角片三个顶点中是否有两个顶点与最小顶点p_Min和最大顶点p_Max重合：如果重合，则该三角片与上一三角片相邻，转向步骤3)，否则转向步骤1)。

3) 判断该三角片的第三个顶点 (不是最小和最大顶点) 的z值与zValue的大小，如果大于zValue，则将该顶点设为最大顶点p_Max，原先的最小顶点p_Min不变；否则，将其设为最小顶点p_Min，原先的最大顶点p_Max不变；更新数据。

4) 通过式 (4) 计算切平面与该三角片的下一个交点，记入交点链表。

5) 标记该三角片并剔除，避免重复循环检索，以提高生成轮廓线数据的速度。

6) 转向步骤1)，继续检索，直到回到第一个相交点p₁，得到一条有向封闭的轮廓线。

3 实验系统设计及结果

针对以上的基本思路，建立如图 5所示的实验系统，首先从计算机层析成像 (Computed Tomography, CT) 或者核磁共振成像 (Magnetic Resonance Imaging, MRI) 获取人体组织图像数据；然后通过医学图像三维重建面绘制算法生成三角片模型，该模型本质上和STL模型是一致的，都是由点、线、面以及法向量组成；最后对生成的三维模型进行切平面分层算法计算分层轮廓线，以实现快速成型加工。

本文共采用三组CT断层图像数据：分别是145幅胸椎，像素分辨率为125×105；133幅脊柱，像素分辨率为100×100；149幅头骨，像素分辨率为256×256。应用Intel Xeon CPU E5-2643 3.30 GHz处理器，32 GB内存，通过VS2010和OpenGL混合编程。不同组织对应的医学图像序列的重建模型如图 6所示。

表 1显示了每个骨组织对应的实验参数，其中三角片数目是指重建后模型的三角片总数，分层时间对应的分层厚度设置为0.1 mm，STL模型分别是由各组织图像数据对应生成的。从表 1可以看出，分组后的分层时间相对于STL模型分层时间有所减少提高，平均效率提高了4.65%。

图 7依次显示了计算后的胸椎、头骨和脊柱分层轮廓。

计算得到三角片模型轮廓线后，通过文献[17]提到的填充算法将其填充，从而得到3D打印的运行轨迹。随机抽取两层头骨轮廓的填充轨迹在Matlab中仿真如图 8所示。

采用极光尔沃Z-603S和聚乳酸 (Polylactice Acid, PLA) 材料分别对上述三种骨组织进行打印，最终的快速成型结果如图 9所示。从左往右依次是胸椎、头骨和脊柱的3D打印原型。

4 结语

医疗快速成型是一个多学科交叉的领域，在医疗诊断、治疗、手术、医学教育、器官植入和制作手术工具等方面扮演着重要角色。针对医疗3D在骨组织的应用关键技术研究，本文提出了一种从医学图像序列生成骨组织模型轮廓线的简化过程方法。该方法结合医学图像三维重建技术和快速成型技术，将医学图像重建后的三角片模型按照建立顺序分组后与切平面求交点轮廓，减少了切平面遍历三角片数量，与医学图像数据生成的STL模型相比，在分层效率上平均提高了4.65%。通过填充轮廓线可以直接生成3D打印轨迹，并通过3D打印机制作图 9所示的实例，证实了该方法的可行性。但是这种方法也有一定的不足之处，与专业的软件相比，每次制作面绘制模型时都要对程序中所要设定的等值面阈值以及针对不同的图像格式作一些调整。下一步研究重点为克服手动调试，实现自动根据图像信息生成三维模型轮廓以及进一步完善实验系统。