Marching Cube算法在点云重建上的简单应用

Marching Cube算法在点云重建上的简单应用

一、 MC算法的简介

  Marching Cubes算法是面显示算法中的经典算法，处理的对一般是断层扫描（CT），或是核磁共振成像（MRI）等产生的图像。MC算法由Lorensen 和Cline 两人在Siggraph Proceedings (pp. 163-169) 提出。

二、MC算法原理

  算法的基本思想是逐个处理数据场中的立方体（体素），分类出与等值面相交的立方体，采用插值计算出等值面与立方体边的交点。根据立方体每一顶点与等值面的相对位置，将等值面与立方体边的交点按一定方式连接生成等值面，作为等值面在该立方体内的一个逼近表示。之所以这样，是由于Marching Cubes有个基本假设：沿六面体边的数据场呈连续性变化。也就是讲，如果一条边的两个顶点分别大于或小于等值面的值，则在该条边上有且仅有一点是这条边与等值面的交点。  具体的来讲，在二维平面中，对于某棱边，如果它的两个端点v1、v2标记不同，那么等值面一定与此棱边相交。且交点坐标为：              **P=P1+(isovalue-V1)(P2-P1)/(V2-V1)**  其中P代表等值点坐标，P1、P2代表两个端点的坐标，V1、V2代表两个端点的权值（或图像的灰度值），isovalue代表阈值。对于每个四边形来说，每个顶点两种情况（大于或小于），4个顶点共16种情况（图1a），考虑到旋转对称性，从新分类后可得4种基本模式（图1b）。

 在三维空间中，由于每一立方体共有8个顶点，每个顶点共有2个状态（物体内和物体外），因此共有256种组合状态，分析立方体体素的2种对称性：  （1）顶点状态反转，等值三角面片的拓扑结构不变，也就是讲，大于等值面与小于等值面的点是可以相互替换的。  （2）旋转对称性，经过适当旋转，有许多状态是一致的。这样，可归纳出15种模式（见图2）。

    在实现时，可按照立方体顶点状态构造等值面连接模式的查找表，并可直接由立方体各顶点的状态检索出其中等值面的分布模式，确定该立方体体素内的等值面三角片连接方式。 

三、MC算法的步骤

① 根据对称关系构建一个256种相交关系的索引表。该表指明等值面与体素的哪条边相交。② 提取立方体的8个顶点，构成一个体素并把这8个顶点编号。   ③ 根据每个顶点与阈值的比较确定该顶点在面内还是面外。④ 把这8个顶点构成的01串组成一个8位的索引值。⑤ 用索引值在上边的索引表里查找对应关系，并求出与立方体每条边的点。⑥ 用交点构成三角形面片或者是多边形面片。⑦ 遍历三维图像的所有体素，重复执行②到⑥。

四、MC在点云重建上的简单应用

  MC算法主要运用在断层扫描（CT），或是核磁共振成像（MRI）等产生的图像中，但是在点云重建中MC算法也可以产生不错的效果。主要思路就是将点云包围在一个大立方体中，然后将立方体均匀分割为若干个小立方体，把小立方体作为体素来进行处理，其中阈值设为0，则建立的等值面即为点云的表面。该方法主要部分有边表（相交关系的索引表）、三角形表（由边表查找生成的三角形等值面）的构建、体素顶点标记、权值以及法向的求取、求边的等值点。下面分别详细的介绍每部分的过程。  1、构建边表和三角形表  本文构造的边表和三角形表并没有经过简化，而是将256种情况都列举了出来。例如下图立方体就是要处理的一个体素，为了方便讲解对其顶点进行编号：

    边表存放的是等值面与立方体相交关系的索引表，上面讲到一共有2^8种情况，则边表大小为256。立方体有12条边，用12位二进制来表示每条边的具体情况。该边标志为1则表示边上有与等值面的交点，为0则没有。然后用八位二进制表示顶点的情况，为1则表示在面内，为0则表示在面外。则由八位二进制数据可以对应到256种情况的边表中。    将顶点编号后对边也进行编号：

边 编号 ⑧⑦ 1 ⑦⑥ 2 ⑥⑤ 3 ⑤⑧ 4 ④③ 5 ③② 6 ②① 7 ④① 8 ⑧④ 9 ⑦③ 10 ⑥② 11 ⑤① 12

     顶点的八位二进制表示中是从**左边**起从顶点1到顶点8，边的12位二进制是从**右边**起从边1到边12。例如0000 0001则表示顶点8在面内，其余顶点则在面外，由此可知等值面与边⑧④、⑤⑧、⑧⑦相交，则对应的12位二进制表示为0001 0000 1001，化为16进制为0x109。以顶点的八位二进制值为边表偏移，以对应的12位二进制值为边表元素，则上述例子可表示为边表偏移为1的元素值为0x109。然后对256种情况一一对应即可构造出边表。按照本文的方法构造的边表如下：

    const int mcubes::edgeTable[256] = {    0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,    0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,    0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,    0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,    0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,    0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,    0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,    0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,    0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,    0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,    0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,    0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,    0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,    0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,    0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,    0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,    0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,    0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,    0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,    0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,    0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,    0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,    0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,    0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,    0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,    0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   };

    三角形表的构造和边表差不多，也是由顶点的八位二进制对应到三角形表中，所以三角形表也有256种情况，但是三角表形的元素存放的是每种对应情况的三角形的顶点。由MC算法易知体素中最多可以形成5个三角形的等值面，所以本文的三角形表的数据结构为:                  int triTable[256][16]；    例如triTable[1]={0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}，即表示顶点情况为0000 0001（偏移的二进制表示）的等值面的三角形情况。其中-1表示此处不是顶点数据，不做处理，0表示边1上的等值点，8表示边9上的等值点，3表示边4上的等值点。由边表可以快速的构建出三角形查找表，当然也可以不构造边表而直接构造三角形表。  2、求体素每个顶点的标记、权值以及法向    求体素顶点的标记、权值以及法向，需要求出以R(参数)为半径，以顶点为球心的球内的所有点云数据。记顶点为P，邻近点为Ki（i从1到n)，邻近点的法向量为Ni(i从1到n)，邻近点的权重Wi(i从1到n)，顶点的法向量为G,d(p,Ki)表示顶点到Ki点的欧式距离,kiP表示向量。则：

Wi=(1−2.5∗(d(P,Ki)));
value=(∑ni=1KiP−→−∗Ni−→∗Wi)/(∑ni=1Wi)
G→=(∑ni=1Ni−→∗Wi)/(∑ni=1Wi)

    其中value即是顶点的权值，顶点P的标记则由value的正负来判定，和点云的法向量的整体朝向有关。  3、求边的等值点           由于Marching Cubes有个基本假设：沿六面体边的数据场呈连续性变化。所以边的等值点由插值可以求出。具体代码如下：

Point mcubes::cacuInterval(Point p1, Point p2) { //计算交点，参数为边的两个顶点;    Point p;    double diff;//权重；    diff = (isovalue - p1.v) / (p2.v - p1.v);//isovalue阈值，P.v是顶点P的权值;    p.x = p1.x + (p2.x - p1.x) * diff;//等值点x坐标;    p.y = p1.y + (p2.y - p1.y) * diff;//等值点y坐标;    p.z = p1.z + (p2.z - p1.z) * diff;//等值点z坐标;    p.nx = p1.nx + (p2.nx - p1.nx) * diff;//等值点法向量x分量;    p.ny = p1.ny + (p2.ny - p1.ny) * diff;//等值点法向量y分量;    p.nz = p1.nz + (p2.nz - p1.nz) * diff;//等值点法向量z分量;    float len = sqrt(p.nx*p.nx+p.ny*p.ny+p.nz*p.nz);//等值点法向量单位化;    p.nx /= len;    p.ny /= len;    p.nz /= len;    return p;}

五、测试结果

  测试数据是bunny数据，txt格式，存储内容为点坐标及法向量。测试结果如下：

     测试结果显示效果不错，当然算法实现简单，在效率上还有待改进。