基于八叉树的网格生成算法剖析

前言

　　对于网格生成这个主题，之前的网格生成系列的三篇博客文章分别介绍了MC算法，SMC算法以及Cuberille算法三种方法。同时还有一篇介绍网格生成与种子点生长算法高效结合的算法。本篇文章继续这一主题，介绍采用八叉树结构来进行网格生成的算法，这个算法并不是独立于之前介绍的算法，而是基于了SMC算法，同时也采纳了一些MC算法范畴内的思想。换句话说，就是使用八叉树对上述网格生成进行改良，从而进一步减少网格的规模。

研究动机与SMC算法的网格规模

　　在SMC算法那一篇中已经指出，SMC算法是相对于MC算法的简化改良，它的一大特点就是生成网格的规模相对于MC算法大大减少。但是SMC算法也毕竟是从体元中产生三角片，三角片的大小是不会超过一个体元的范围的。而三维图像往往会有很大的尺寸这就意味着体元的数量会非常大。因而即便不采用MC算法而采用SMC算法，仍然会产生大量密集的小三角片。比如下图是Engine数据的三维图像模型，可以看出三角网格是由很多小三角片构成的。



　　进一步观察这样的模型，我们发现：模型中实际上存在大量的共面三角片。了解网格型的话就会意识到这不是一个好的网格模型，因为网格模型通过小三角形片去拟合模型表面的曲面。而像上图那样的平面，完全可以使用更少的三角片去表示。比如一个具有下图左方结构的Mesh用8三角片个表示一个正方形，完全等价于右图用2个三角片表示的正方形。

8个三角片组成一个正方形	2个也能组个一样的

　　因此这就引入了进一步减少网格规模的主题，这也是本文使用八叉树来简化网格生成的研究动机。无论是SMC算法还是MC算法，其生成的网格都有存在大量共面三角形的可能，而这些数量众多的小三角形也有合并成较少的大三角形的可能。尤其是SMC算法，由于其小三角形的法向和形状具有良好的可枚举性，因而在局部非常容易出现能够合并的小三角形。

　　为此还必须指出，本文即将讨论的方法与针对Mesh进行网格削减的一系列算法（如之前的文章提到的顶点聚簇）并不是同一类算法。后者旨在脱离图像的概念，仅从几何意义上去对Mesh进行削减。而本文的算法是在网格生成之前，通过八叉树结构进行网格生成的优化，从而直接在Mesh被构造时就减少网格的规模。所以其本质上还是算作Mesh Generation算法而不是Mesh Processing算法。

空间划分树的分类-均等树和BON树

　　熟悉计算机的都知道树这个数据结构，尤其是二叉树。在计算机科学中，树形结构经常被用来组织和检索数据。当然树的用处不止是这些。在计算机图形学中，经常使用树来对空间进行划分，也就是把空间组织成一个层次结构。这样树的每一个节点就被赋予了空间上的意义，可以用来代表一个区域，而树的父子关系也正好能表示区域之间的包含关系。一般的，一维空间可以使用二叉树来划分，而平面和空间区域分别使用四叉树和八叉树来进行划分。

　　对于一维的数据，比如一个下面的范围[0-12]，可以使用二叉树划分。最终使得范围内的每一个子元素都在二叉树的节点上。这样的划分方式也比较常见，简单的对范围空间除以2即可。所以扩展到平面和空间就是对各个轴坐标范围不停的除以2，这种方式与快速排序中递归树的划分方式有一点相似（注意不是完全一样的），我们把这种方式划分的树叫做均等树。比如下图展示了对0-12范围的数据进行均等划分到底所建立的树。均等树对于[A,B]这样的范围，首先找到(A+B)/2，然后就将区域分成[A,(A+B)/2]，[(A+B)/2,B]。之后再继续按相同的方式划分到不能再分为止。

均等树的划分方式

　　BON树是一种不大同于上述均等树划分方式的树。从上图可以看出均等树的算法在除2的时候会遇到奇数偶数的问题。实际上在遇到奇数个元素的范围内，均等树所划分的子树实际上是不均等的，按上文说的办法，遇到奇数个元素，就会把中位数分在靠前的范围中，而靠后的范围实际上会少一个元素。这样划分到最后会出现单枝叶子节点的情况。而BON树为了保证树形，会在一开始就把最初的范围变成2的幂，从而能让划分的范围一直处于偶数个元素。例如下图展示使用BON树的分法来划分0-12的范围，首先需要找到包含这个范围的最小的2的幂，即16，然后再对这个范围进行划分，超出范围的节点不再创建。

BON树实际能表示到最近的2的幂	BON树超出部分不建立分支

　　BON树的划分逻辑相比于均等树，并不更加复杂，而且相比于均等树存在一个潜在的好处。例如上图中对0-12范围所创建的均等树，可以从节点的编号上就能够获得节点的位置信息。例如9这个节点，9的二进制位为1001，正好对应这这个9的节点被BON树所分的层次信息。从图上0对应左子树，1对应右子树，这样从根节点到1001这个数就正好是右→左→左→右。这样同时提供了树节点插入的一个思路：在表示0-M范围的BON树中，若要插入A，则可以跟据A的二进制位和树的层数来确定这个A在树中的位置。同时，BON树的非叶子节点都能指代一个范围。例如9上面的父节点，就能代表100X（X为0或1）这个由两个数组成的范围，而这个父节点的父亲又能指代10XX的范围。



　　以上是对两种具有不同空间划分方式的树的介绍，本文所要继续介绍的八叉树划分方法，是基于BON树的。同时在下文中会详细说明这样的树结构究竟是如何与具体的网格生成算法结合的。

树结构与信息的概括

　　在介绍具体算法之前，还需要继续对空间划分树的一大特点进行介绍，而对这个特点的利用也是本文算法的基本思想。这个特点就是树结构具有对局部信息的概括能力。

　　了解算法的人都比较清楚，像树这种具有层次的数据结构，很好的实现了对信息的概括。比如二叉排序树组织数据，相比与使用数组来说，树的每一个节点都概括了自己子树的信息，就是：左子树下面的都比自己小，右子树下面的都比自己大，这就是一种概括。那么当有数据来检索的时候，与一个节点的比较，就能判断出数组访问一个元素更多的大小关系。这就是层次结构之所以能实现高效检索的原因。而数组的每一个位置不像树那样具有概括信息的能力，访问每个位置都只能知道这一个位置的元素的信息。

　　而对于本文说的空间划分树，可以通过下面的例子来说明其是如何概括信息的。例如下图中展示了一个编号为0-12的范围，这个范围分布有红绿两种颜色。也就是有的编号处是红色，有的是绿色。



　　使用空间划分树来对这个范围组织一个BON树结果如下：



　　假如我们想用最少的节点来表示这样一个红绿数组，那么就可以对这个树执行一个shrink操作。Shrink代表收缩是expand的反义词，类似于我们在treeview菜单上HEAD的那个节点上点击了”-”的动作一样。

Expand后的树	Shrink后的树

　　这样每一个节点都会概括自己的子节点，当发现他们的颜色一致的时候，就会把这些信息概括到自己身上。下图展示了这种自底向上概括的过程，最后形成的树充分概括了这段红绿数组的信息。

	原始BON树
	shrink倒数第一层
	shrink倒数第二层
	shrink倒数第三层，不能再shrink完毕

　　当对最后shrink完毕的树检索到最右的节点时，在这个节点上能够得知他所表示的范围即8-12都是绿色。

　　了解了树的自底向上的shrink操作，就会在大方向上清楚了本文的算法将如何减少三角形的面片数。如下图所示，如果左图中这样八个体元中的三角形都被抽取出来，这样形成的一个大三角形平面是由4个三角形组成的，而其实如果能够将这些体元合并后再进行抽取。那么抽取出的三角形就只一个。这样实现了用更少的三角形表示相同的平面区域。这也就是本文所述的算法所要达到的直接目的。

相邻的八个体元中的三角片	体元合并在一起后再抽取，只有一个三角片了

基于八叉树的网格生成

　　在介绍MC算法和SMC算法的时候我们已经知道，三角形片都是从边界体元中抽取出来的。三维图像中一般都有大量的实体元和空体元，这些体元对三角形抽取是无用的，只有边界体元是我们需要的。MC和SMC算法执行过程中，对三维图像的三个轴的三重循环扫描过程实际上就是一个搜寻边界体元的过程，找出所有体元配置不为0和255的体元，就是边界体元，每找到这样的体元，就抽取其中的三角形。在种子点生长算法与网格生成算法结合的那文章里面，也是这样的方式，无非就是寻找边界体元的方式变成基于种子点生长的方式而不是全部扫描。所以我们可以把使用MC、SMC算法生成网格的步骤总结成三步：

寻找边界体元；
抽取体元三角片；
组合三角片成Mesh；

　　而本文的基于八叉树的网格算法，通过树的shrink来减少需要抽取三角片的体元数量，这样在上面步骤的基础上有如下的改变：

寻找边界体元；
将边界体元插入八叉树；
对八叉树执行shrink；
抽取shrink后的各体元三角片；
组合三角片成Mesh；

　　从上文的步骤可以想到，shrink之前树中存着都是同一层的等大小的体元，而shrink之后，就一部分节点被父节点吸收概括成一个超体元，那么这棵树里的体元就会由各种大小不同的体元构成，既有边长为1的单位体元，也有边长为2的幂如2、4、8等的超体元。

　　上文的步骤是粗略的一个步骤，其中最为关键的是中间三步，每一步都还有具体需要关心的问题：第二步将边界体元插入八叉树，是如何的插法；第三步对树进行shrink，那么符合什么样条件的节点可以shrink，shrink到什么程度为止；第四步中提取体元的三角片，如果是单位体元，提取的方式和SMC算法一样，但如果是超体元将如何提取三角片。解决好了上述问题，就完整实现了本文所介绍的算法。

八叉树的创建以及边界体元的插入

　　首先需要解决建树以及插入节点的问题。本算法的实现使用的是BON八叉树，空间树的建立必然关联着一个空间范围，边界体元集合是有范围的，而且至多也不会超出三维图像的范围。为了建树，首先定义树的节点类型OctreeNode，声明其为模版类型方便于携带不同的参数，更具一般性。

public class OctreeNode<T>
{
public OctreeNode<T>[] Children;//孩子指针,数组大小为8
public OctreeNode<T> Parent;//父节点指针
public T Parms;//携带的参数
public int XMin;//所代表范围的X轴下界
public int YMin;//所代表范围的Y轴下界
public int ZMin;//所代表范围的Z轴下界
public int XMax;//所代表范围的X轴下界
public int YMax;//所代表范围的Y轴下界
public int ZMax;//所代表范围的Z轴下界
public int IndexInParent;//自己在父节点孩子数组中的索引
public int LayerIndex;//自己所在的层索引
public bool IsLeaf()
{
return (XMin == XMax)&&(YMin==YMax)&&(ZMin==ZMax);
}//返回是否是叶子节点
public OctreeNode()
{
}
public override string ToString()
{
if (IsLeaf())
{
return string.Format("[{0},{1}][{2},{3}][{4},{5}] {6} Leaf", XMin, XMax, YMin, YMax, ZMin, ZMax,Parms);
}
if (Parms == null)
return string.Format("[{0},{1}][{2},{3}][{4},{5}] {6}", XMin, XMax, YMin, YMax, ZMin, ZMax, "not simple"); ;
return string.Format("[{0},{1}][{2},{3}][{4},{5}] {6}", XMin, XMax, YMin, YMax, ZMin, ZMax,Parms);
}
}//BON八叉树节点

　　上述定义中为了方便附加了很多信息如节点的各轴范围，节点的层次以及在父节点的索引等，这些信息有的是可以动态获取的，开辟空间记录下来是用空间换时间，本文为方便实现，附加了比较多的这些辅助信息。

　　那么下一步就是定义一颗BON树RegionOctree：

public class RegionOctree<T>
{
private static int GetMax2Power(int xmax,int ymax,int zmax,ref int log)
{
int max = xmax;
if (ymax > max)
max = ymax;
if (zmax > max)
max = zmax;
if ((max & (max - 1)) == 0)
{
double L = Math.Log(max, 2);
log = (int)L + 1;
return max;
}
else
{
double L = Math.Log(max, 2);
log = (int)L + 2;
return (int)Math.Pow(2, log - 1);
}
}
private int Width;//树所关联空间范围的X上界
private int Height;//树所关联空间范围的Y上界
private int Depth;//树所关联空间范围的Z上界
public OctreeNode<T> Root;//树根节点
public  int NodeCount;//所有节点总数
public  int LeafCount;//叶子节点
private int Scale;//2的幂包围盒边长
private int LayerNum;//层次数
private OctreeNode<T>[] NodeLayers;//指代一条由根通往叶子的路径

public RegionOctree(int width,int height,int depth)//使用范围构造BON树
{
this.Width = width;
this.Height = height;
this.Depth = depth;
Scale = GetMax2Power(Width,Height,Depth,ref LayerNum);
NodeCount = 0;
Root = new OctreeNode<T>();
Root.XMin = 0;
Root.XMax = Scale-1;
Root.YMin = 0;
Root.YMax = Scale-1;
Root.ZMin = 0;
Root.ZMax = Scale-1;
Root.Parent = null;
Root.IndexInParent = -1;
Root.LayerIndex = LayerNum - 1;
Root.Children = new OctreeNode<T>[8];
NodeLayers = new OctreeNode<T>[LayerNum];
NodeLayers[0] = Root;
}

public OctreeNode<T> CreateToLeafNode(int x,int y,int z)
{
LeafCount++;
for (int i = 1; i <= LayerNum - 1; i++)
{
int index = GetIndexOn(x, y, z, LayerNum - i-1);
if (NodeLayers[i - 1].Children[index] == null)
{
OctreeNode<T> node = new OctreeNode<T>();
NodeCount++;
node.Parent = NodeLayers[i - 1];
node.IndexInParent = index;
node.Children = new OctreeNode<T>[8];
node.LayerIndex = NodeLayers[i - 1].LayerIndex - 1;
InitRangeByParentAndIndex(node, NodeLayers[i - 1], index);
NodeLayers[i - 1].Children[index] = node;
}
NodeLayers[i]=NodeLayers[i-1].Children[index];
}
return NodeLayers[NodeLayers.Length - 1];
}//将关联着坐标（x，y，z）处元素一路插入到底层为叶子节点

private int GetIndexOn(int x, int y, int z, int bitindex)
{
int ret = 0;
if ((x & (1 << bitindex)) != 0)
{
ret |= 1;
}
if ((y & (1 << bitindex)) != 0)
{
ret |= 2;
}
if ((z & (1 << bitindex)) != 0)
{
ret |= 4;
}
return ret;
}

private void InitRangeByParentAndIndex(OctreeNode<T> node,OctreeNode<T> pnode, int index)
{
int deltaX = (pnode.XMax - pnode.XMin + 1) / 2;
int deltaY = (pnode.YMax - pnode.YMin + 1) / 2;
int deltaZ = (pnode.ZMax -  pnode.ZMin + 1) / 2;
if ((index & 1) == 0)
{
node.XMin = pnode.XMin;
node.XMax = pnode.XMin + deltaX - 1;
}
else
{
node.XMin = pnode.XMin + deltaX;
node.XMax = pnode.XMax;
}
if ((index & 2) == 0)
{
node.YMin = pnode.YMin;
node.YMax = pnode.YMin + deltaY - 1;
}
else
{
node.YMin = pnode.YMin + deltaY;
node.YMax = pnode.YMax;
}
if ((index & 4) == 0)
{
node.ZMin = pnode.ZMin;
node.ZMax = pnode.ZMin + deltaZ - 1;
}
else
{
node.ZMin = pnode.ZMin + deltaZ;
node.ZMax = pnode.ZMax;
}
}//使用父节点的信息初始化子节点的范围
}

　　从上述代码可以看出，对三维图像范围(0-width,0-height,0-depth)关联一颗BON树，那么这颗BON树就首先需要明确自己的层数，也就是需要找到最小的2的幂的包围盒。GetMax2Power函数实现了这一功能，例如一个具有范围(200,300,400)的空间，寻找到的最小2的幂包围盒就是(512,512,512)。

　　在理解插入方法之前还需要对子节点的顺序进行固定的编号，在这里采用的方式是按0-7的二进制位000-111进行分配：从右数第一位1表示X轴正方向；0表示X轴负方向。从右数第二位1表示Y轴正方向；0表示Y轴负方向。从右数第三位1表示Z轴正方向；0表示Z轴负方向。所以8个子体元的编号位置对应关系就如下表所示。

	X范围：（XMIN~XMAX） Y范围：（YMIN~YMAX） Z范围：（ZMIN~ZMAX） 设各方向范围中点分别为：XMID,YMID,ZMID	孩子编号 对应范围 0 （XMIN-XMID，YMIN-YMID，ZMIN-ZMID） 1 （XMID-XMAX，YMIN-YMID，ZMIN-ZMID） 2 （XMIN-XMID，YMID-YMAX，ZMIN-ZMID） 3 （XMID-XMAX，YMID-YMAX，ZMIN-ZMID） 4 （XMIN-XMID，YMIN-YMID，ZMID-ZMAX） 5 （XMID-XMAX，YMIN-YMID，ZMID-ZMAX） 6 （XMIN-XMID，YMID-YMAX，ZMID-ZMAX） 7 （XMID-XMAX，YMID-YMAX，ZMID-ZMAX）
图示	父节点信息	孩子信息

　　为了理解与插入相关的三个函数CreateToLeafNode、GetIndexOn与InitRangeByParentAndIndex。下面使用一组数据来举例说明一个具有坐标(80,85,22)的元素是如何插入一个范围为(200,200,210)的BON树的。首先根据这三个轴坐标的最大值210，找到对应比它大的最小的2的幂为256，所以这个BON树所能表示的的最大范围是(256，256，256)而层数（或者说深度）为8（256为2的8次方）。

　　确定了BON树的层数为8后，将这(80,85,22)的三个坐标的二进制位像下图一样填在3个长度为8（等于层数）的表里，我们就可以从这个表中得出插入这个节点的方法：



　　图中的层内索引是每一竖列的二进制位组合起来的结果，也是X行的值*4与Y行的值*2与Z行的值*1相加的结果。这个层内索引表达了这个节点在每一层应当被插入到哪一颗子树。从上图的结果看，层内索引分别为第0、第3、第0、第7、第0、第6、第2。这就意味着将(80,85,22)插入这课BON树的步骤为：

将[80,85,22]插入0层第0个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入1层第3个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入2层第0个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入3层第7个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入4层第0个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入5层第6个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入6层第4个孩子节点，若该孩子节点为NULL则创建之。
将[80,85,22]插入7层第2个孩子节点，若该孩子节点为NULL则创建之。

　　CreateToLeafNode函数主逻辑便是这样一个过程，这样的过程结束后，私有成员NodeLayers内的信息包含了所有经过的节点。GetIndexOn函树通过检查位来找到对应层的层内索引。InitRangeByParentAndIndex函数则负责初始化每一个节点的创建信息。例如一个父节点所表示的范围为[X0-X1,Y0-Y1,Z0-Z1]，那么他的子节点所表示的范围就正好将每个轴的范围平均劈成两半。对应的信息可以参见上文的表中孩子信息中的范围。

　　这样八叉树中如何插入边界体元信息的问题就解决了，体元集合的插入就相当于是三维坐标的依次插入，这样插入之后每一个边界体元都处在八叉树的叶子节点上。当然这个八叉树只是可能的最大分支数为8。一般来说，会有很多节点不满8个子节点。所以这样的树就叫做Branch On Need树，简称为BON树。

Shrink操作的合并原则

　　在之前介绍算法目标的时候就已经提到，本算法的目的是减少共面三角形数，也就是企图合并那些在一个平面上的小三角片，将其合成大三角片。如下图所示的4个例子情况，都是可以合并的情况：

合并前
合并后

　　从上图中可见，假如一个节点的子节点里面的三角片全部是共平面的，这样的节点就可以shrink。再进一步考虑，有的体元配置中只有一个平面，而有的体元却包含了了多个不共面的三角片，如下图所示：

不能shrink的例子之一	不能shrink的例子之二	能shrink的例子之一	能shrink的例子之二

　　那么后者就绝对没有被shrink的可能。而前者的体元若作为某节点的孩子，假如他的兄弟节点和他一样具有共面的三角片，那么这些兄弟就可以与之合并成一个大的超体元，那么现在的问题就是如何判断三角形的共面情况。考虑到SMC算法的三角形片形状和方向是可以枚举的那几种形状，我们可以把这些信息结合体元配置做一个映射，从而快速判断三角形片是否是共面的。使用空间解析几何的知识我们可以知道：空间中的所有平面都可以使用方程Ax+By+Cz=D来表示。其中（A,B,C）的组合表征平面的方向，而D值表征位移，确定了这四个系数，平面就被唯一确定。

　　在SMC算法介绍中，我们已经知道三角形片的顶点都是体元的体素。下面以单个体元为例子来说明三角片的平面方程。下图是笔者使用WPF技术制作的指示SMC体元配置和三角片关系的软件截图，可以通过http://files.cnblogs.com/chnhideyoshi/Debug.rar下载。这个软件通过拖动滑块来观看不同的体元配置的三角片情况。

共面配置之一
共面配置之二
共面配置之三
非共面配置

　　我们可以看出，凡是具有多个不共面三角形的情况，都被标注为“非共面配置”；而只有一个正三角形或者是由两个三角形组合成长方形这样的一类的体元配置被标注为“共面配置”。具有共面配置的体元就存在着与兄弟节点合并的可能。同时还可以发现，这些三角形的方程的形式是有限的，ABC的值只在0、1、-1三个数中变化。D值的变化也是有限的，但由于这只是基准点坐标为（0,0,0）的体元，一旦将这个体元坐标一般化为（Dx，Dy，Dz），那么这个D值就会跟随这个坐标进行变化，而ABC是不变的。

　　所谓共面配置只是描述了体元内三角形的形式，或者说三角形平面方向，即方程中的(ABC)三元组，而三角形所在平面的具体位置还依赖于D，D不同即使ABC相同，这样的三角形也只能是平行的三角形，而不是共面的。所以能够合并的体元，其中的三角形不仅要具有相同的法向，也要有相同的D值。也就是说，八叉树的节点可以记录下所代表体元中三角形片的法向和D值，这两个信息关系到体元能否与兄弟节点合并。

　　所以综合以上的陈述，总结一下对于八叉树的一个节点P，他的孩子C0-C7能够shrink到P必须遵循以下的原则：

C0-C7要么为NULL（即非边界体元区域），要么为一种共面配置。
C0-C7不为NULL的对应体元的共面配置必须均为同一种
C0-C7的不为NULL的节点均有相同的D值
C0-C7的后代节点若不为NULL则需均不含合并失败的节点，也就是一个节点若shrink失败，其父节点不再shrink。

　　由于使用三元组来表达三角片方向这个信息会有诸多冗余，考虑到法向量是有限可枚举的，所以可以为这些法向量进行编号，设置查找数组，以编号来代替三元组来进行存储。设置的查找数组如下：

public struct Int32Quad
{
public int A;
public int B;
public int C;
public int D;
public Int32Quad(int a, int b, int c, int d)
{
this.A = a;
this.B = b;
this.C = c;
this.D = d;
}
public override string ToString()
{
return string.Format("{0}x+{1}y+{2}z={3}", A, B, C, D);
}
}//代表一个平面方程对象
public class OctreeTable
{
public static byte[] ConfigToNormalTypeId = new byte[256]
{
13,0,1,2,3,13,4,13,5,6,13,13,7,13,13,8,3,9,13,13,13,13,13,13,13,13,13,0,13,13,13,13,5,13,10,13,13,13
,13,1,13,13,13,13,13,13,13,13,7,13,13,11,13,13,13,13,13,13,13,13,13,13,13,2,0,13,13,13,9,13,13,13,13
,13,13,13,13,13,3,13,13,13,13,13,13,13,13,13,13,13,
13,13,13,13,13,13,6,13,13,13,13,13,12,13,13,13,13,13,13,13,13,4,13,13,5,13,13,13,13,10,13,13,13,13
,13,13,13,1,1,13,13,13,13,13,13,13,10,13,13,13,13,5,13,
13,4,13,13,13,13,13,13,13,13,12,13,13,13,13,13,6,13,13,13,13,13,13,
13,13,13,13,13,13,13,13,13,13,13,3,13,13,13,13,13,13,13,13,13,9,13,13,13,0,2,13,13,13,
13,13,13,13,13,13,13,13,11,13,13,7,13,13,13,13,13,13,13,13,1,13,13,13,13,10,13,5,13,13,
13,13,0,13,13,13,13,13,13,13,13,13,9,3,8,13,13,7,13,13,6,5,13,4,13,3,2,1,0,13
};//体元配置对应的法向量索引
public static Int16Triple[] NormalTypeIdToNormal = new Int16Triple[13]
{
new Int16Triple(1,-1,-1),
new Int16Triple(1,-1,1),
new Int16Triple(1,-1,0),
new Int16Triple(1,1,1),
new Int16Triple(1,0,1),
new Int16Triple(1,1,-1),
new Int16Triple(1,0,-1),
new Int16Triple(1,1,0),
new Int16Triple(1,0,0),
new Int16Triple(0,1,1),
new Int16Triple(0,1,-1),
new Int16Triple(0,1,0),
new Int16Triple(0,0,1)
};//枚举的法向量集合

public static byte[] ConfigToEqType = new byte[256]
{
55,0,1,2,3,55,4,55,5,6,55,55,7,55,
55,8,9,10,55,55,55,55,55,55,55,
55,55,11,55,55,55,55,12,55,13,55,
55,55,55,14,55,55,55,55,55,55,55
,55,15,55,55,16,55,55,55,55,55,55
,55,55,55,55,55,17,18,55,55,55,19,
55,55,55,55,55,55,55,55,55,20,55,
55,55,55,55,55,55,55,55,55,55,55,
55,55,55,55,0,21,55,55,55,55,55,22,
55,55,55,55,55,55,55,55,23,55,55,24,55,
55,55,55,25,55,55,55,0,55,0,0,26,27,55,55,
55,55,55,55,55,28,55,55,55,55,29,55,55,
30,55,55,55,55,55,55,55,55,31,55,55,55,
55,55,32,55,55,55,55,55,55,55,55,55,55,
55,55,55,55,55,0,55,33,55,55,55,55,
55,0,55,55,55,34,55,0,0,35,36,55,55,55,
55,55,55,55,55,55,55,55,37,55,55,38,55,
55,55,55,55,55,55,0,39,55,55,0,55,40,0,41,55,
55,55,55,42,55,55,0,55,55,55,0,55,0,43,44,45,55
,55,46,55,0,47,48,55,49,0,50,51,52,53,55
};//体元配置对应的方程索引，55表示非共面配置
public static Int32Quad[] EqTypeToEqQuad = new Int32Quad[54]
{
new Int32Quad(1,-1,-1,-1),
new Int32Quad(1,-1,1,0),
new Int32Quad(1,-1,0,0),
new Int32Quad(1,1,1,1),
new Int32Quad(1,0,1,1),
new Int32Quad(1,1,-1,0),
new Int32Quad(1,0,-1,0),
new Int32Quad(1,1,0,1),
new Int32Quad(1,0,0,1),
new Int32Quad(1,1,1,2),
new Int32Quad(0,1,1,1),
new Int32Quad(1,-1,-1,0),
new Int32Quad(1,1,-1,1),
new Int32Quad(0,1,-1,0),
new Int32Quad(1,-1,1,1),
new Int32Quad(1,1,0,1),
new Int32Quad(0,1,0,0),
new Int32Quad(1,-1,0,1),
new Int32Quad(1,-1,-1,0),
new Int32Quad(0,1,1,1),
new Int32Quad(1,1,1,2),
new Int32Quad(1,0,-1,0),
new Int32Quad(0,0,1,1),
new Int32Quad(1,0,1,2),
new Int32Quad(1,1,-1,0),
new Int32Quad(0,1,-1,-1),
new Int32Quad(1,-1,1,2),
new Int32Quad(1,-1,1,1),
new Int32Quad(0,1,-1,0),
new Int32Quad(1,1,-1,1),
new Int32Quad(1,0,1,1),
new Int32Quad(0,0,1,0),
new Int32Quad(1,0,-1,1),
new Int32Quad(1,1,1,1),
new Int32Quad(0,1,1,0),
new Int32Quad(1,-1,-1,1),
new Int32Quad(1,-1,0,0),
new Int32Quad(0,1,0,1),
new Int32Quad(1,1,0,2),
new Int32Quad(1,-1,1,0),
new Int32Quad(0,1,-1,1),
new Int32Quad(1,1,-1,2),
new Int32Quad(1,-1,-1,-1),
new Int32Quad(0,1,1,2),
new Int32Quad(1,1,1,3),
new Int32Quad(1,0,0,0),
new Int32Quad(1,1,0,0),
new Int32Quad(1,0,-1,-1),
new Int32Quad(1,1,-1,-1),
new Int32Quad(1,0,1,0),
new Int32Quad(1,1,1,0),
new Int32Quad(1,-1,0,-1),
new Int32Quad(1,-1,1,-1),
new Int32Quad(1,-1,-1,-2),
};//体元配置对应的三角片方程集合
public const byte NormalNotSimple = 13;
}

　　这样我们就可以为八叉树节点Node的T类型成员指定一个类型NodeParms，这个类型的变量携带体元的共面信息NormalTypeId和D，为方便起见还加上了体元配置Config，一共占据2个字节加一个整型的空间。以便在shrink执行时提取该信息进行判断。

public class NodeParms
{
public byte Config;
public byte NormalTypeId;
public int D;
public override string ToString()
{
if (NormalTypeId < OctreeTable.NormalTypeIdToNormal.Length)
return "\"" + NormalTypeId + "," + D + "\"";
else
return "not simple";
}
}

　　这里附上完整的OctreeSurfaceGenerator的代码如下：

public static class OctreeSurfaceGenerator
{
#region consts
public static byte VULF = 1 << 0;
public static byte VULB = 1 << 1;
public static byte VLLB = 1 << 2;
public static byte VLLF = 1 << 3;
public static byte VURF = 1 << 4;
public static byte VURB = 1 << 5;
public static byte VLRB = 1 << 6;
public static byte VLRF = 1 << 7;
//以上为体素为实点的位标记
public static Int16Triple[] PointIndexToPointDelta = new Int16Triple[8]
{
new Int16Triple(0, 1, 1 ),
new Int16Triple(0, 1, 0 ),
new Int16Triple(0, 0, 0 ),
new Int16Triple(0, 0, 1 ),
new Int16Triple(1, 1, 1 ),
new Int16Triple(1, 1, 0 ),
new Int16Triple(1, 0, 0 ),
new Int16Triple(1, 0, 1 )
};//体元内每个体素相对基准体素坐标的偏移
public static byte[] PointIndexToFlag = new byte[8]
{
VULF,
VULB,
VLLB,
VLLF,
VURF,
VURB,
VLRB,
VLRF
};//每个体素对应的位标记
#endregion
public static Mesh GenerateSurface(BitMap3d bmp)
{
int width = bmp.width;
int height = bmp.height;
int depth = bmp.depth;
Int16Triple[] tempArray = new Int16Triple[8];
#region CreateTree
RegionOctree<NodeParms> otree = new RegionOctree<NodeParms>(width, height, depth);
Queue<OctreeNode<NodeParms>> nodequeue = new Queue<OctreeNode<NodeParms>>();
for (int k = 0; k < depth-1; k++)
{
for (int j = 0; j < height-1; j++)
{
for (int i = 0; i < width-1; i++)
{
byte value = 0;
for (int pi = 0; pi < 8; pi++)
{
tempArray[pi].X = i + PointIndexToPointDelta[pi].X;
tempArray[pi].Y = j + PointIndexToPointDelta[pi].Y;
tempArray[pi].Z = k + PointIndexToPointDelta[pi].Z;
if (InRange(bmp, tempArray[pi].X, tempArray[pi].Y, tempArray[pi].Z)&&
IsWhite(bmp, tempArray[pi].X, tempArray[pi].Y, tempArray[pi].Z))
{
value |= PointIndexToFlag[pi];
}
}
if (value != 0 && value != 255)
{
OctreeNode<NodeParms> leafnode = otree.CreateToLeafNode(i, j, k);
leafnode.Parms = new NodeParms();
leafnode.Parms.Config = value;
leafnode.Parms.NormalTypeId = OctreeTable.ConfigToNormalTypeId[value];
leafnode.Parms.D = CaculateDFromNormalAndCoord(i, j, k, value);
nodequeue.Enqueue(leafnode.Parent);
}
}
}
}
#endregion
#region Shrink
while (nodequeue.Count != 0)
{
OctreeNode<NodeParms> node = nodequeue.Dequeue();
byte normalType = OctreeTable.NormalNotSimple;
int D = int.MinValue;
if (CanMergeNode(node, ref normalType, ref D))
{
node.Parms = new NodeParms();
node.Parms.NormalTypeId = normalType;
node.Parms.D = D;
nodequeue.Enqueue(node.Parent);
}
}
#endregion
#region ExtractTriangles
MeshBuilder_IntegerVertex mb = new MeshBuilder_IntegerVertex(width+1, height+1, depth+1);
nodequeue.Enqueue(otree.Root);
while (nodequeue.Count != 0)
{
OctreeNode<NodeParms> node = nodequeue.Dequeue();
if (node.Parms == null)
{
for (int i = 0; i < 8; i++)
{
if (node.Children[i] != null)
nodequeue.Enqueue(node.Children[i]);
}
}
else
{
if (node.Parms.NormalTypeId != OctreeTable.NormalNotSimple)
{
if (node.IsLeaf())
{
GenerateFaceLeaf(node, mb, ref tempArray, bmp);
}
else
{
GenerateFace(node, mb, ref tempArray, bmp);
}
}
else
{
if (node.IsLeaf())
{
GenerateFaceLeaf(node, mb, ref tempArray, bmp);
}
else
{
for (int i = 0; i < 8; i++)
{
if (node.Children[i] != null)
nodequeue.Enqueue(node.Children[i]);
}
}
}
}
}//采用层次遍历寻找需要抽取三角片的节点
#endregion
return mb.GetMesh();
}

private static bool InRange(BitMap3d bmp, int x, int y, int z)
{
return x > 0 && x < bmp.width && y > 0 && y < bmp.height && z > 0 && z < bmp.depth;
}
private static bool IsWhite(BitMap3d bmp, int x, int y, int z)
{
return bmp.GetPixel(x, y, z) == BitMap3d.WHITE;
}

private static void GenerateFace(OctreeNode<NodeParms> node, MeshBuilder_IntegerVertex mb, ref Int16Triple[] tempArray, BitMap3d bmp)
{
InitVoxelPositionForNodeRange(node.XMin, node.XMax, node.YMin, node.YMax, node.ZMin, node.ZMax, ref tempArray);
//需要找到该节点的端点位置
byte cubeConfig = 0;
for (int pi = 0; pi < 8; pi++)
{
if (InRange(bmp, tempArray[pi].X, tempArray[pi].Y, tempArray[pi].Z)
&& IsWhite(bmp,tempArray[pi].X,tempArray[pi].Y ,tempArray[pi].Z))
{
cubeConfig |= PointIndexToFlag[pi];
}
}
int index = 0;
while (MCTable.TriTable[cubeConfig, index] != -1)
{
int ei1 = MCTable.TriTable[cubeConfig, index];
int ei2 = MCTable.TriTable[cubeConfig, index + 1];
int ei3 = MCTable.TriTable[cubeConfig, index + 2];

Int16Triple p1 = GetIntersetedPointAtEdge(node, ei1, OctreeTable.NormalTypeIdToNormal[node.Parms.NormalTypeId], node.Parms.D);
Int16Triple p2 = GetIntersetedPointAtEdge(node, ei2, OctreeTable.NormalTypeIdToNormal[node.Parms.NormalTypeId], node.Parms.D);
Int16Triple p3 = GetIntersetedPointAtEdge(node, ei3, OctreeTable.NormalTypeIdToNormal[node.Parms.NormalTypeId], node.Parms.D);

mb.AddTriangle(p1, p2, p3);
index += 3;
}
}//对非叶子节点的超体元的抽取需要参考MCTable求取被截断边的信息

private static void GenerateFaceLeaf(OctreeNode<NodeParms> node, MeshBuilder_IntegerVertex mb, ref Int16Triple[] tempArray, BitMap3d bmp)
{
for (int k = 0; k < 8; k++)
{
tempArray[k].X = node.XMin + PointIndexToPointDelta[k].X;
tempArray[k].Y = node.YMin + PointIndexToPointDelta[k].Y;
tempArray[k].Z = node.ZMin + PointIndexToPointDelta[k].Z;
}
byte value = node.Parms.Config;
int index = 0;
while (SMCTable.TableFat[value, index] != -1)
{
Int16Triple t0 = tempArray[SMCTable.TableFat[value, index]];
Int16Triple t1 = tempArray[SMCTable.TableFat[value, index + 1]];
Int16Triple t2 = tempArray[SMCTable.TableFat[value, index + 2]];
mb.AddTriangle(t0, t1,t2);
index += 3;
}
}//对叶子节点的单位体元的抽取和SMC算法中的抽取一致

private static bool CanMergeNode(OctreeNode<NodeParms> node, ref byte normalType, ref int D)
{
for (int i = 0; i < 8; i++)
{
if (node.Children[i] != null)
{
if (node.Children[i].Parms == null)
{
return false;//说明其下存在不能合并的节点 合并失败
}
else
{
if (node.Children[i].Parms.NormalTypeId != OctreeTable.NormalNotSimple)
{
normalType = node.Children[i].Parms.NormalTypeId;
D = node.Children[i].Parms.D;//记录其中的共面配置信息
}
else
{
return false;//遇到了非共面配置 合并失败
}
}
}
}
for (int i = 0; i < 8; i++)
{
if (node.Children[i] != null)
{
if (node.Children[i].Parms.NormalTypeId != normalType || node.Children[i].Parms.D != D)
{
return false;//体元配置均为共面类型但不是同一种的话也不中
}
}
}
return true;
}

private static Int16Triple GetIntersetedPointAtEdge(OctreeNode<NodeParms> node, int edgeIndex, Int16Triple normal, int d)
{
int x = 0, y = 0, z = 0;
switch (edgeIndex)
{
case 0: { x = node.XMin; y = node.YMax + 1; return new Int16Triple(x, y, (d - normal.X * x - normal.Y * y) / normal.Z); }
case 2: { x = node.XMin; y = node.YMin; return new Int16Triple(x, y, (d - normal.X * x - normal.Y * y) / normal.Z); }
case 4: { x = node.XMax + 1; y = node.YMax + 1; return new Int16Triple(x, y, (d - normal.X * x - normal.Y * y) / normal.Z); }
case 6: { x = node.XMax + 1; y = node.YMin; return new Int16Triple(x, y, (d - normal.X * x - normal.Y * y) / normal.Z); }

case 8: { y = node.YMax + 1; z = node.ZMax + 1; return new Int16Triple((d - normal.Y * y - normal.Z * z) / normal.X, y, z); }
case 9: { y = node.YMax + 1; z = node.ZMin; return new Int16Triple((d - normal.Y * y - normal.Z * z) / normal.X, y, z); }
case 10: { y = node.YMin; z = node.ZMin; return new Int16Triple((d - normal.Y * y - normal.Z * z) / normal.X, y, z); }
case 11: { y = node.YMin; z = node.ZMax + 1; return new Int16Triple((d - normal.Y * y - normal.Z * z) / normal.X, y, z); }

case 1: { x = node.XMin; z = node.ZMin; return new Int16Triple(x, (d - normal.X * x - normal.Z * z) / normal.Y, z); }
case 3: { x = node.XMin; z = node.ZMax + 1; return new Int16Triple(x, (d - normal.X * x - normal.Z * z) / normal.Y, z); }
case 5: { x = node.XMax + 1; z = node.ZMin; return new Int16Triple(x, (d - normal.X * x - normal.Z * z) / normal.Y, z); }
case 7: { x = node.XMax + 1; z = node.ZMax + 1; return new Int16Triple(x, (d - normal.X * x - normal.Z * z) / normal.Y, z); }

default: throw new Exception();
}
}

private static int CaculateDFromNormalAndCoord(int cx, int cy, int cz, byte config)
{
byte index = OctreeTable.ConfigToEqType[config];
if (index >= OctreeTable.EqTypeToEqQuad.Length)
return int.MinValue;
Int32Quad eq = OctreeTable.EqTypeToEqQuad[index];
return eq.D + eq.A * cx + eq.B * cy + eq.C * cz;
}

private static void InitVoxelPositionForNodeRange(int xmin, int xmax, int ymin, int ymax, int zmin, int zmax, ref Int16Triple[] temp)
{
temp[0].X = xmin;  //(0, 1, 1, VULF);
temp[0].Y = ymax + 1;
temp[0].Z = zmax + 1;

temp[1].X = xmin;//(0, 1, 0, VULB);
temp[1].Y = ymax + 1;
temp[1].Z = zmin;

temp[2].X = xmin;//(0, 0, 0, VLLB);
temp[2].Y = ymin;
temp[2].Z = zmin;

temp[3].X = xmin;//(0, 0, 1, VLLF);
temp[3].Y = ymin;
temp[3].Z = zmax + 1;

temp[4].X = xmax + 1;//(1, 1, 1, VURF);
temp[4].Y = ymax + 1;
temp[4].Z = zmax + 1;

temp[5].X = xmax + 1; //(1, 1, 0, VURB);
temp[5].Y = ymax + 1;
temp[5].Z = zmin;

temp[6].X = xmax + 1;//(1, 0, 0, VLRB);
temp[6].Y = ymin;
temp[6].Z = zmin;

temp[7].X = xmax + 1;  //(1, 0, 1, VLRF);
temp[7].Y = ymin;
temp[7].Z = zmax + 1;
}

public static void Test1()
{
BitMap3d image = BitMap3d.CreateSampleTedVolume(100);
Mesh m=GenerateSurface(image);
PlyManager.Output(m, "D://VTKproj//Tree_Test6_T1.ply");
}

public static void Test2()
{
BitMap3d image = BitMap3d.CreateSampleEngineVolume();
Mesh m = GenerateSurface(image);
PlyManager.Output(m, "D://VTKproj//Tree_Test6_T2.ply");
}

public static void Test()
{
Test2();
}
}

　　可以看出在主函数的函数体中的三大步逻辑：

创建树
Shrink
三角形提取

　　在函数CanMergeNode中充分体现出了本段所述的合并原则。

体元三角形的抽取

　　虽然已经贴出了完整的代码，但是最后还剩下一个需要解答的问题，就是如何抽取体元三角形。首先在SMC算法中我们知道如果是单位体元，直接根据体元配置查找三角形表，这样就可以找到三角形的顶点和连接方式。三角形的顶点一定是体元的8个体素点之一。而在本文介绍的算法中，涉及到了对非单位体元、也就是超体元的三角片提取。不难想到，八叉树的叶子节点存储的都是单位体元的坐标，在shrink之后若其没有被上一级吸收则仍然为单位体元，那么这样的节点抽取三角形仍然采用SMC算法的办法；但是对于超体元中三角片的提取，则有必要描述一下这之中的技巧。

　　对于超体元首先必须明确这样一个事实：虽然超体元中的三角形片的顶点也一定是体素点，但这些顶点不见得是超体元的8个端点，而是可能是其下层节点体元的端点。如下图所示的节点为叶子节点的父亲，包含了8个单位体元。可以看出，合并成一个超体元后，三角片穿过这个超体元的边，而不一定经过这个超体元的端点。
　　所以这里在提取三角形的时候，就不再是只考虑使用SMC三角形表，因为SMC三角形表只能提取单位体元的三角形。由于MC三角形表涉及到了对三角形穿过体元边的信息，所以这里对超体元提取三角形还得用上MC的三角形表。

　　在提取三角形之前还需要考虑超体元配置的求法，由于超体元本质上表示了一个立方体范围的体元坐标，所以其中包含了多个体元的体素，我们需要考察超体元端点的体素来决定这个超体元的配置。下面的表通过对比来显示超体元和单位体元求取体元配置的不同

体素编号 体素名称 体素位置 0 VULF （x，y+1，z+1） 1 VULB （x，y+1，z） 2 VLLB （x，y，z） 3 VLLF （x，y，z+1） 4 VURF （x+1，y+1，z+1） 5 VURB （x+1，y+1，z） 6 VLRB （x+1，y，z） 7 VLRF （x+1，y，z+1）	体素编号 体素名称 体素位置 0 VULF （xmin，ymax+1，zmax+1） 1 VULB （xmin，ymax+1，zmin） 2 VLLB （xmin，ymin，zmin） 3 VLLF （xmin，ymin，zmax+1） 4 VURF （xmax+1，ymax+1，zmax+1） 5 VURB （xmax+1，ymax+1，zmin） 6 VLRB （xmax+1，ymin，zmin） 7 VLRF （xmax+1，ymin，zmax+1）
单位体元，坐标为（x,y,z）	超体元，指代范围（xmin-xmax,ymin-ymax,zmin-zmax）

　　在求取超体元的插值点的时候，根据边的方向和三角形的方程，可以枚举出每条边上插值点如下表所示：

边代号 边两端顶点索引 边方程 与三角片方程AX+BY+CZ=D的交点位置 e0 （0，1） x=xmin,y=ymax+1 （xmin,ymax+1,(D-A*xmin-B*(ymax+1))/C） e1 （1，2） x=xmin,z=zmin （xmin,(D-A*xmin-C*zmin)/B,zmin） e2 （2，3） x=xmin,y=ymin （xmin,ymin,(D-A*xmin-B*ymin)/C） e3 （3，0） x=xmin,z=zmax+1 （xmin,(D-A*xmin-C*(zmax+1))/B,zmax+1） e4 （4，5） x=xmax+1,y=ymax+1 （xmax+1,ymax+1,(D-A*(xmax+1)-B*(ymax+1))/C） e5 （5，6） x=xmax+1,z=zmin （xmax+1,(D-A*(xmax+1)-C*zmin)/B,zmin） e6 （6，7） x=xmax+1,y=ymin （xmax+1,ymin,(D-A*(xmax+1)-B*ymin)/C） e7 （7，4） xmax+1,z=zmax+1 （xmax+1,(D-A*(xmax+1)-C*(zmax+1))/B,zmax+1） e8 （0，4） y=ymax+1,z=zmax+1 （(D-B*(ymax+1)-C*(zmax+1))/A,ymax+1,zmax+1） e9 （1，5） y=ymax+1,z=zmin （(D-B*(ymax+1)-C*zmin)/A,ymax+1,zmin） e10 （2，6） y=ymin,z=zmin （(D-B*ymin-C*zmin)/A,ymin,zmin） e11 （3，7） y=ymin,z=zmax+1 （(D-B*ymin-C*(zmax+1))/A,ymin,zmax+1）

图示[b]超体元，指代范围（xmin-xmax,ymin-ymax,zmin-zmax）[/b]

边信息

　　这样就解释了第三个问题，同时也说明了RegionOctree类的其他函数的含义。

算法实现与测试

　　除了上文贴出的代码之外，还涉及到哈希表类、Mesh类等在之前的博文中都有涉及到，这里不再重复贴出。

　　算法结果的测试采用两组数据，一组是TED三维图像，其中的内容部分表现为一个三棱锥。还有一组是Engine数据。算法测试时会与SMC算法的结果进行对比。

　　首先是算法的局部结果对比：

TED数据	SMC结果预览
八叉树算法结果预览

Engine数据	SMC结果预览
八叉树算法结果预览

　　算法结果参数对比如下表所示：

	SMC算法顶点数	八叉树算法顶点数	SMC算法面片数	八叉树算法面片数
TED数据	16471	3123	32938	5147
Engine数据	216147	144117	432370	259479

　　可见八叉树算法针对有大片平面区域的模型，能够一定程度上减少网格的规模。算法工程的下载地址为：https://github.com/chnhideyoshi/SeededGrow2d