膝关节镜手术是一种微创手术，具有创伤小、恢复快等优点.但是由于膝关节腔空间狭小，内窥镜视野范围小，手术难度高^[1].得益于医学成像技术的发展，个性化患者数据的采集比以往更容易，数据的精确度更高.如何利用这些数据应用于定量化的术前规划和术中导航，是一个值得关注的课题.

目前对于腔体提取的研究大多是针对规则、密闭内腔的三维模型，如周敬之等^[2]利用特征识别重构管道模型空腔；胡庆夕和吴福等^[3]利用射线法检测内腔边界实现内腔的提取，但这不适用于开放的膝关节腔.根据膝关节解剖结构的特点，本文提出一种基于区域生长法的自动提取膝关节腔体的方法，为进一步对腔体环境进行量化分析提供基础.

图 1所示的是由CT图像三维重建得到的膝关节的股骨和胫骨模型^[4]，它们是围合膝关节腔体的两个最重要的结构.因此，膝关节腔体自动提取的一个重要目标就是识别并提取股骨的下端面和胫骨的上端面.

图 1(Figure 1)

图 1 膝关节三维模型 Figure 1 The three-dimensional model of the knee

本文所提出的腔体自动提取算法，其流程可以概括为：首先，分别在股骨和胫骨上选取一个属于腔体边界面的三角面片作为种子，使用区域生长法初步提取两个边界面；然后对提取到的两个边界面作孔洞修补和曲面边缘光顺处理；最后对腔体的上下两个边界面进行缝合处理.

种子面的识别，需要利用法线与相邻解剖结构的三角面片的相交关系判断.例如图 2中，胫骨上的一个三角面片法线n₁的外法线方向与股骨的三角面片相交，则该三角面片属于腔体边界面.因此，在股骨(胫骨)上选取种子面片，要计算三角面片的法线是否与胫骨(股骨)上的三角面片相交.同时，在区域生长过程中也要将这一条件作为约束.

图 2(Figure 2)

图 2 腔体边界面的法线相交情况 Figure 2 Intersection of the normal and the triangular facets

由于膝关节模型的面片数量大，若对所有的三角面片逐个测试，计算量非常大.为了减小计算量，本文引入了轴向包围盒(AABB), 并建立相应的八叉树数据结构.进行相交测试时，本文依据八叉树逐层遍历，这样就避免了对每个三角面片逐个相交测试，提高计算效率.

一个给定对象的轴向包围盒被定义为能包含该对象且各边平行于坐标轴的最小正六面体^[5].如图 3所示，将包围盒分为8个子六面体，对子六面体依次编号为0~7.然后对子六面体的三维空间区域逐个卦限测试，如果子六面体包含的面片数量大于预先设定的值则继续分解，直至子六面体包含的面片数小于设定值.假设其中的6号子六面体所包含的面片数依然大于设定值，则将6号子六面体分解.

图 3(Figure 3)

图 3 轴对齐包围盒(AABB) Figure 3 Axis aligned bounding box

在计算法线是否与空间三角面片相交时，本文使用Möller和Trumbore提出的方法，将向量表示成一个基点和方向向量^[6]：

可以简单地将一个三角形定义为一组有序顶点{V₀, V₁, V₂}，如图 4所示.

图 4(Figure 4)

图 4 直线与空间三角形的相交 Figure 4 Intersection of the line and spatial triangle

三角形内的任何点都可以用它相对于三角形的顶点的位置来定义：

其中u+v+w=1.三元组(u, v, w)称为Q的重心坐标.由于w=1-(u+v)，因此这里仅使用(u, v)来表示重心坐标.将方程(1) 代入(2) 中，就能计算法线——三角形的相交：$P + t\mathit{\boldsymbol{\hat d}} = \left( {1-\left( {u + v} \right)} \right){V_0} + u{V_1} + v{V_2}$, 该式可以扩展为

这是一个由具有3个未知数的方程所组成的线性系统，可以用克莱姆法则^[7]求解：

一旦求得了t、u和v，就能通过考察这3个参数的值来确定交点是否位于三角形内(而不是位于三角形所在平面的其他地方)：如果0≤u≤1，0≤v≤1并且u+v≤1，那么相交于三角形内；否则就不相交于三角形内.

区域生长法首先从原始面片数据中选取一个初始三角面片，并作为种子面加入三角曲面集，3条边分别加入边界集^[8].然后按照一定的规则从边界处向边界曲面集添加新的三角面片，以此向外生长成一个完整的曲面.

图 5描述了区域增长法的原理：首先选取一个满足前述条件的三角面片作为种子；然后以种子面片(黑色填充面)为起点，分别测试种子面片的3个邻接面(蓝色填充面)是否满足腔体面片的条件；将满足条件的面片添加到目标三角面片集合；接着以新添加的面片为种子，迭代向外生长^[9].

图 5(Figure 5)

图 5 区域生长法生长过程 Figure 5 The processing of region growing

图 6是根据区域生长提取出的两个曲面，由于外法线相交于邻接结构的三角面片只是腔体边界面片的充分条件，并不是所有的腔体边界面片都满足这一条件，所以以上述条件为约束并不能将所有的腔体边界面片提取出来，就形成了图 7中的一些孔洞.

图 6(Figure 6)

图 6 腔体边界面的初步提取效果 Figure 6 Preliminary extraction effect of the boundary surface

图 7(Figure 7)

图 7 腔体边界面上的孔洞 Figure 7 The holes of the boundary surface

判别孔洞的依据是曲面内部的边界.如图 8所示，如果三角曲面中的一条边只属于一个三角面片，则称之为边界边^[10-11].把所有邻接的边界边顺次连接所形成的封闭环即为孔洞边界(内环，位于曲面的内部)或者是曲面的外边界(外环)^[12-13].

图 8(Figure 8)

图 8 孔洞边界信息 Figure 8 Boundary of the hole

假设曲面的法线方向为穿过纸面向外，并设每条边都是有方向的，根据右手定则，每个三角形的环向都是逆时针方向(如图 8右图所示).由一序列有向边构成的外环也是逆时针方向的，而内环却是顺时针方向的，以此便可以区分内环和外环，从而可以判断某个环是曲面的外边界还是孔洞的边界.

若检测出某环线是曲面的内环，则采用区域生长法以环线的边界三角形为种子，遍历原始三角面片数据，若检测到三角面片与该边界三角形邻接，便将其添加至目标三角面片集合，以此迭代向孔洞内生长即可修补孔洞.

图 9为孔洞修补后的效果图.通过上述方法提取的腔体曲面边呈锯齿状，需要对它进行光顺处理.

图 9(Figure 9)

图 9 孔洞修补后的效果 Figure 9 The effect of filling holes

边缘光顺处理的方法是迭代地对满足条件的的凸角(见图 10)进行删除，或者对凹角进行填补，直到曲面边缘环线上的最小凸(凹)角的角度大于本文设定的阈值.

图 10(Figure 10)

图 10 外环上的凸(凹)角 Figure 10 The convex and concave angle of outer ring

在边缘轮廓中，相邻两条边的夹角越小，意味着轮廓线在该处的变化很大.因此，在迭代的每一步，都需要沿着外环环线(曲面的边缘轮廓)方向遍历，计算每个夹角的余弦值，找出余弦值最大的夹角(即角度最小的角)；然后以右手定则判断该角为凸角还是凹角.若该角为凹角，则连接凹角顶点的前一顶点和后一顶点，增加一个三角形，将凹角补齐；若该角为凸角，则将该顶点连接的三角形删去.直至遍历环线上的最小凸(凹)角角度大于本文设定的阈值(最小角度的余弦值).经过实验验证，将阈值设为-0.55既可以达到很好的光顺效果，又能保持其原有的轮廓.图 11为光顺处理后的效果.

图 11(Figure 11)

图 11 曲面边缘光顺后的效果 Figure 11 The effect of edge smoothing

从股骨(上)和胫骨(下)上提取出两个腔体边界面并作光顺处理后，还要将这两曲面缝合成封闭的曲面.

由于股骨和胫骨之间的空间很小，提取出的两曲面之间的距离很小，所以在缝合时只需要添加包含上下曲面边缘轮廓顶点的三角面即可达到缝合的效果.如图 12所示，设折线V₁V₂V₃是上表面的外环线上的一段，P₁P₂P₃是下表面的外环线上的一段.首先找出环1上到P₁距离最小的顶点；然后连接V₁P₁和V₁P₂构成一个三角面V₁P₁P₂.以同样的方法可以增加V₂P₂P₃和V₂P₃P₄和其他连接环线1上的一个顶点和环线2上的两个顶点的三角面；最后，只要反向连接环线1上的两个顶点和环线2上的一个顶点，即可补全两表面的侧面，得到封闭的腔体曲面如图 13所示.

图 12(Figure 12)

图 12 侧面三角面片的增加 Figure 12 Suturing surfaces by adding triangles

图 13(Figure 13)

图 13 曲面缝合后的效果 Figure 13 Suturing of the surfaces

本文通过分析膝关节腔的解剖结构的特点，确定判别腔体边界面片的一个充分条件：待测三角面的外法线与邻近结构的三角面相交.利用该条件可以确定种子面和判断待测三角面是否属于边界三角面.通过区域生长法自动提取出腔体的初步轮廓曲面，并对这些轮廓面进行孔洞修复、边缘光滑和曲面缝合处理.实验证明，该算法可以自动地提取并重建膝关节腔体模型.此外，针对算法中耗时最大的环节——外法线与三角面的相交测试，利用八叉树这种存储结构，大大地减少了法线与面片的相交计算，从而提高了计算的效率.