流体模拟（三）

Marching Cube算法（1）

我们在实现流体表面重建时，需要事先在空间中划分网格，我们的流体系统正好已经完成了此项工作。其次利用Marching Cube算法计算出构成表面的三角片面，最后根据所有三角片面的计算每个顶点的法线值。

Marching Cube算法，最早是Lorensen等人在“Marching cubes: A high resolution 3D surface construction algorithm”这篇文章中提出的算法，他们最早是用来实现医学图像的三维重建，现在该算法可以用于任何基于标量场的表面重建，并还有许多优化版本。本文就提一下它的最初版本。

我们给空间划分为了许多小网格，每个网格也称为体元（cell），如下图所示：

由于该算法是基于标量场的，所以我们需要知道网格点是在场外还是场内（表面外还是表面内），即f(x,y,z)=0的点为表面点，所以也称点(x,y,z)满足该隐函数。由于每个体元上有八个顶点。若连续的两个顶点，有一个顶点位于表面外，另一个顶点位于表面内，则我们的表面就肯定位于这两个顶点之间，即可以用加权插值获得最后表面的位置。因此我们的算法需要遍历所有的网格，判断网格的每个顶点是在表面外还是表面内，若为表面内则赋值1，否则赋值0。

Lorensen等人他们就发现了每个网格一共就8个顶点，那所有的可能性就只有 $2^{8}$ =256种。并且由于立方体具有旋转对称性，经过一定的旋转会得到同样的值，因此他们经过归纳发现所有情况可以归结为以下15种：

之后给每个立方体定义了一个固定的顶点编号和边编号：

如图所示，v为顶点编号e为边编号。因为正好是八个顶点，就可以利用一个char类型（8位）来保存该编号索引。由于我们最后的三角片面的顶点都会在边上，立方体有12条边，所以我们可以用一个12位2进制数表示边是否含有三角片面的点。因为我们构建的三角形片面情况一共只有256种，所以我们可以建立边索引表。如我们顶点索引位00000001，说明只有v1在表面内，因此说明表面的顶点在e1，e4，e9上，则边索引为001000001001，化为16进制为0x109。以顶点的八位二进制值为边表索引，以对应的12位二进制值为边表值，则上述例子可表示为边表索引为1的值为0x109。然后对256种情况一一对应即可构造出边表。构造的边表如下：

const unsigned int 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
};

构造完边表之后，我们还需要构造一个三角形表。三角形表记录了体元中三角形片面的连接方式。根据那15种情况，可以发现一个体元内的三角面最多只有4个，所以我们的三角形表可以定义为一个256*12的2维数组。同样通过顶点索引出三角形连接方式。如同样的0000001。则triTable[1]={0, 8, 3, X, X, X, X, X, X, X, X, X, X, X, X, X}（X定义为255即不存在点）。表示当三角形在顶点v1在表面内时，只有一个三角形，且三角形的三个点在边e1，e9，e4上。我们贴上三角形表的代码：

#define X 255
const int triTable[256][16] = {{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X},{0, 8, 3, X, X, X, X, X, X, X, X, X, X, X, X, X},{0, 1, 9, X, X, X, X, X, X, X, X, X, X, X, X, X},{1, 8, 3, 9, 8, 1, X, X, X, X, X, X, X, X, X, X},{1, 2, 10, X, X, X, X, X, X, X, X, X, X, X, X, X},{0, 8, 3, 1, 2, 10, X, X, X, X, X, X, X, X, X, X},{9, 2, 10, 0, 2, 9, X, X, X, X, X, X, X, X, X, X},{2, 8, 3, 2, 10, 8, 10, 9, 8, X, X, X, X, X, X, X},{3, 11, 2, X, X, X, X, X, X, X, X, X, X, X, X, X},{0, 11, 2, 8, 11, 0, X, X, X, X, X, X, X, X, X, X},{1, 9, 0, 2, 3, 11, X, X, X, X, X, X, X, X, X, X},{1, 11, 2, 1, 9, 11, 9, 8, 11, X, X, X, X, X, X, X},{3, 10, 1, 11, 10, 3, X, X, X, X, X, X, X, X, X, X},{0, 10, 1, 0, 8, 10, 8, 11, 10, X, X, X, X, X, X, X},{3, 9, 0, 3, 11, 9, 11, 10, 9, X, X, X, X, X, X, X},{9, 8, 10, 10, 8, 11, X, X, X, X, X, X, X, X, X, X},{4, 7, 8, X, X, X, X, X, X, X, X, X, X, X, X, X},{4, 3, 0, 7, 3, 4, X, X, X, X, X, X, X, X, X, X},{0, 1, 9, 8, 4, 7, X, X, X, X, X, X, X, X, X, X},{4, 1, 9, 4, 7, 1, 7, 3, 1, X, X, X, X, X, X, X},{1, 2, 10, 8, 4, 7, X, X, X, X, X, X, X, X, X, X},{3, 4, 7, 3, 0, 4, 1, 2, 10, X, X, X, X, X, X, X},{9, 2, 10, 9, 0, 2, 8, 4, 7, X, X, X, X, X, X, X},{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, X, X, X, X},{8, 4, 7, 3, 11, 2, X, X, X, X, X, X, X, X, X, X},{11, 4, 7, 11, 2, 4, 2, 0, 4, X, X, X, X, X, X, X},{9, 0, 1, 8, 4, 7, 2, 3, 11, X, X, X, X, X, X, X},{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, X, X, X, X},{3, 10, 1, 3, 11, 10, 7, 8, 4, X, X, X, X, X, X, X},{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, X, X, X, X},{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, X, X, X, X},{4, 7, 11, 4, 11, 9, 9, 11, 10, X, X, X, X, X, X, X},{9, 5, 4, X, X, X, X, X, X, X, X, X, X, X, X, X},{9, 5, 4, 0, 8, 3, X, X, X, X, X, X, X, X, X, X},{0, 5, 4, 1, 5, 0, X, X, X, X, X, X, X, X, X, X},{8, 5, 4, 8, 3, 5, 3, 1, 5, X, X, X, X, X, X, X},{1, 2, 10, 9, 5, 4, X, X, X, X, X, X, X, X, X, X},{3, 0, 8, 1, 2, 10, 4, 9, 5, X, X, X, X, X, X, X},{5, 2, 10, 5, 4, 2, 4, 0, 2, X, X, X, X, X, X, X},{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, X, X, X, X},{9, 5, 4, 2, 3, 11, X, X, X, X, X, X, X, X, X, X},{0, 11, 2, 0, 8, 11, 4, 9, 5, X, X, X, X, X, X, X},{0, 5, 4, 0, 1, 5, 2, 3, 11, X, X, X, X, X, X, X},{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, X, X, X, X},{10, 3, 11, 10, 1, 3, 9, 5, 4, X, X, X, X, X, X, X},{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, X, X, X, X},{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, X, X, X, X},{5, 4, 8, 5, 8, 10, 10, 8, 11, X, X, X, X, X, X, X},{9, 7, 8, 5, 7, 9, X, X, X, X, X, X, X, X, X, X},{9, 3, 0, 9, 5, 3, 5, 7, 3, X, X, X, X, X, X, X},{0, 7, 8, 0, 1, 7, 1, 5, 7, X, X, X, X, X, X, X},{1, 5, 3, 3, 5, 7, X, X, X, X, X, X, X, X, X, X},{9, 7, 8, 9, 5, 7, 10, 1, 2, X, X, X, X, X, X, X},{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, X, X, X, X},{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, X, X, X, X},{2, 10, 5, 2, 5, 3, 3, 5, 7, X, X, X, X, X, X, X},{7, 9, 5, 7, 8, 9, 3, 11, 2, X, X, X, X, X, X, X},{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, X, X, X, X},{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, X, X, X, X},{11, 2, 1, 11, 1, 7, 7, 1, 5, X, X, X, X, X, X, X},{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, X, X, X, X},{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, X},{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, X},{11, 10, 5, 7, 11, 5, X, X, X, X, X, X, X, X, X, X},{10, 6, 5, X, X, X, X, X, X, X, X, X, X, X, X, X},{0, 8, 3, 5, 10, 6, X, X, X, X, X, X, X, X, X, X},{9, 0, 1, 5, 10, 6, X, X, X, X, X, X, X, X, X, X},{1, 8, 3, 1, 9, 8, 5, 10, 6, X, X, X, X, X, X, X},{1, 6, 5, 2, 6, 1, X, X, X, X, X, X, X, X, X, X},{1, 6, 5, 1, 2, 6, 3, 0, 8, X, X, X, X, X, X, X},{9, 6, 5, 9, 0, 6, 0, 2, 6, X, X, X, X, X, X, X},{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, X, X, X, X},{2, 3, 11, 10, 6, 5, X, X, X, X, X, X, X, X, X, X},{11, 0, 8, 11, 2, 0, 10, 6, 5, X, X, X, X, X, X, X},{0, 1, 9, 2, 3, 11, 5, 10, 6, X, X, X, X, X, X, X},{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, X, X, X, X},{6, 3, 11, 6, 5, 3, 5, 1, 3, X, X, X, X, X, X, X},{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, X, X, X, X},{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, X, X, X, X},{6, 5, 9, 6, 9, 11, 11, 9, 8, X, X, X, X, X, X, X},{5, 10, 6, 4, 7, 8, X, X, X, X, X, X, X, X, X, X},{4, 3, 0, 4, 7, 3, 6, 5, 10, X, X, X, X, X, X, X},{1, 9, 0, 5, 10, 6, 8, 4, 7, X, X, X, X, X, X, X},{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, X, X, X, X},{6, 1, 2, 6, 5, 1, 4, 7, 8, X, X, X, X, X, X, X},{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, X, X, X, X},{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, X, X, X, X},{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, X},{3, 11, 2, 7, 8, 4, 10, 6, 5, X, X, X, X, X, X, X},{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, X, X, X, X},{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, X, X, X, X},{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, X},{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, X, X, X, X},{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, X},{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, X},{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, X, X, X, X},{10, 4, 9, 6, 4, 10, X, X, X, X, X, X, X, X, X, X},{4, 10, 6, 4, 9, 10, 0, 8, 3, X, X, X, X, X, X, X},{10, 0, 1, 10, 6, 0, 6, 4, 0, X, X, X, X, X, X, X},{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, X, X, X, X},{1, 4, 9, 1, 2, 4, 2, 6, 4, X, X, X, X, X, X, X},{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, X, X, X, X},{0, 2, 4, 4, 2, 6, X, X, X, X, X, X, X, X, X, X},{8, 3, 2, 8, 2, 4, 4, 2, 6, X, X, X, X, X, X, X},{10, 4, 9, 10, 6, 4, 11, 2, 3, X, X, X, X, X, X, X},{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, X, X, X, X},{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, X, X, X, X},{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, X},{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, X, X, X, X},{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, X},{3, 11, 6, 3, 6, 0, 0, 6, 4, X, X, X, X, X, X, X},{6, 4, 8, 11, 6, 8, X, X, X, X, X, X, X, X, X, X},{7, 10, 6, 7, 8, 10, 8, 9, 10, X, X, X, X, X, X, X},{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, X, X, X, X},{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, X, X, X, X},{10, 6, 7, 10, 7, 1, 1, 7, 3, X, X, X, X, X, X, X},{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, X, X, X, X},{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, X},{7, 8, 0, 7, 0, 6, 6, 0, 2, X, X, X, X, X, X, X},{7, 3, 2, 6, 7, 2, X, X, X, X, X, X, X, X, X, X},{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, X, X, X, X},{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, X},{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, X},{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, X, X, X, X},{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, X},{0, 9, 1, 11, 6, 7, X, X, X, X, X, X, X, X, X, X},{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, X, X, X, X},{7, 11, 6, X, X, X, X, X, X, X, X, X, X, X, X, X},{7, 6, 11, X, X, X, X, X, X, X, X, X, X, X, X, X},{3, 0, 8, 11, 7, 6, X, X, X, X, X, X, X, X, X, X},{0, 1, 9, 11, 7, 6, X, X, X, X, X, X, X, X, X, X},{8, 1, 9, 8, 3, 1, 11, 7, 6, X, X, X, X, X, X, X},{10, 1, 2, 6, 11, 7, X, X, X, X, X, X, X, X, X, X},{1, 2, 10, 3, 0, 8, 6, 11, 7, X, X, X, X, X, X, X},{2, 9, 0, 2, 10, 9, 6, 11, 7, X, X, X, X, X, X, X},{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, X, X, X, X},{7, 2, 3, 6, 2, 7, X, X, X, X, X, X, X, X, X, X},{7, 0, 8, 7, 6, 0, 6, 2, 0, X, X, X, X, X, X, X},{2, 7, 6, 2, 3, 7, 0, 1, 9, X, X, X, X, X, X, X},{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, X, X, X, X},{10, 7, 6, 10, 1, 7, 1, 3, 7, X, X, X, X, X, X, X},{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, X, X, X, X},{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, X, X, X, X},{7, 6, 10, 7, 10, 8, 8, 10, 9, X, X, X, X, X, X, X},{6, 8, 4, 11, 8, 6, X, X, X, X, X, X, X, X, X, X},{3, 6, 11, 3, 0, 6, 0, 4, 6, X, X, X, X, X, X, X},{8, 6, 11, 8, 4, 6, 9, 0, 1, X, X, X, X, X, X, X},{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, X, X, X, X},{6, 8, 4, 6, 11, 8, 2, 10, 1, X, X, X, X, X, X, X},{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, X, X, X, X},{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, X, X, X, X},{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, X},{8, 2, 3, 8, 4, 2, 4, 6, 2, X, X, X, X, X, X, X},{0, 4, 2, 4, 6, 2, X, X, X, X, X, X, X, X, X, X},{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, X, X, X, X},{1, 9, 4, 1, 4, 2, 2, 4, 6, X, X, X, X, X, X, X},{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, X, X, X, X},{10, 1, 0, 10, 0, 6, 6, 0, 4, X, X, X, X, X, X, X},{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, X},{10, 9, 4, 6, 10, 4, X, X, X, X, X, X, X, X, X, X},{4, 9, 5, 7, 6, 11, X, X, X, X, X, X, X, X, X, X},{0, 8, 3, 4, 9, 5, 11, 7, 6, X, X, X, X, X, X, X},{5, 0, 1, 5, 4, 0, 7, 6, 11, X, X, X, X, X, X, X},{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, X, X, X, X},{9, 5, 4, 10, 1, 2, 7, 6, 11, X, X, X, X, X, X, X},{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, X, X, X, X},{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, X, X, X, X},{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, X},{7, 2, 3, 7, 6, 2, 5, 4, 9, X, X, X, X, X, X, X},{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, X, X, X, X},{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, X, X, X, X},{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, X},{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, X, X, X, X},{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, X},{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, X},{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, X, X, X, X},{6, 9, 5, 6, 11, 9, 11, 8, 9, X, X, X, X, X, X, X},{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, X, X, X, X},{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, X, X, X, X},{6, 11, 3, 6, 3, 5, 5, 3, 1, X, X, X, X, X, X, X},{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, X, X, X, X},{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, X},{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, X},{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, X, X, X, X},{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, X, X, X, X},{9, 5, 6, 9, 6, 0, 0, 6, 2, X, X, X, X, X, X, X},{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, X},{1, 5, 6, 2, 1, 6, X, X, X, X, X, X, X, X, X, X},{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, X},{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, X, X, X, X},{0, 3, 8, 5, 6, 10, X, X, X, X, X, X, X, X, X, X},{10, 5, 6, X, X, X, X, X, X, X, X, X, X, X, X, X},{11, 5, 10, 7, 5, 11, X, X, X, X, X, X, X, X, X, X},{11, 5, 10, 11, 7, 5, 8, 3, 0, X, X, X, X, X, X, X},{5, 11, 7, 5, 10, 11, 1, 9, 0, X, X, X, X, X, X, X},{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, X, X, X, X},{11, 1, 2, 11, 7, 1, 7, 5, 1, X, X, X, X, X, X, X},{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, X, X, X, X},{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, X, X, X, X},{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, X},{2, 5, 10, 2, 3, 5, 3, 7, 5, X, X, X, X, X, X, X},{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, X, X, X, X},{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, X, X, X, X},{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, X},{1, 3, 5, 3, 7, 5, X, X, X, X, X, X, X, X, X, X},{0, 8, 7, 0, 7, 1, 1, 7, 5, X, X, X, X, X, X, X},{9, 0, 3, 9, 3, 5, 5, 3, 7, X, X, X, X, X, X, X},{9, 8, 7, 5, 9, 7, X, X, X, X, X, X, X, X, X, X},{5, 8, 4, 5, 10, 8, 10, 11, 8, X, X, X, X, X, X, X},{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, X, X, X, X},{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, X, X, X, X},{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, X},{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, X, X, X, X},{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, X},{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, X},{9, 4, 5, 2, 11, 3, X, X, X, X, X, X, X, X, X, X},{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, X, X, X, X},{5, 10, 2, 5, 2, 4, 4, 2, 0, X, X, X, X, X, X, X},{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, X},{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, X, X, X, X},{8, 4, 5, 8, 5, 3, 3, 5, 1, X, X, X, X, X, X, X},{0, 4, 5, 1, 0, 5, X, X, X, X, X, X, X, X, X, X},{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, X, X, X, X},{9, 4, 5, X, X, X, X, X, X, X, X, X, X, X, X, X},{4, 11, 7, 4, 9, 11, 9, 10, 11, X, X, X, X, X, X, X},{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, X, X, X, X},{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, X, X, X, X},{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, X},{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, X, X, X, X},{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, X},{11, 7, 4, 11, 4, 2, 2, 4, 0, X, X, X, X, X, X, X},{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, X, X, X, X},{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, X, X, X, X},{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, X},{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, X},{1, 10, 2, 8, 7, 4, X, X, X, X, X, X, X, X, X, X},{4, 9, 1, 4, 1, 7, 7, 1, 3, X, X, X, X, X, X, X},{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, X, X, X, X},{4, 0, 3, 7, 4, 3, X, X, X, X, X, X, X, X, X, X},{4, 8, 7, X, X, X, X, X, X, X, X, X, X, X, X, X},{9, 10, 8, 10, 11, 8, X, X, X, X, X, X, X, X, X, X},{3, 0, 9, 3, 9, 11, 11, 9, 10, X, X, X, X, X, X, X},{0, 1, 10, 0, 10, 8, 8, 10, 11, X, X, X, X, X, X, X},{3, 1, 10, 11, 3, 10, X, X, X, X, X, X, X, X, X, X},{1, 2, 11, 1, 11, 9, 9, 11, 8, X, X, X, X, X, X, X},{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, X, X, X, X},{0, 2, 11, 8, 0, 11, X, X, X, X, X, X, X, X, X, X},{3, 2, 11, X, X, X, X, X, X, X, X, X, X, X, X, X},{2, 3, 8, 2, 8, 10, 10, 8, 9, X, X, X, X, X, X, X},{9, 10, 2, 0, 9, 2, X, X, X, X, X, X, X, X, X, X},{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, X, X, X, X},{1, 10, 2, X, X, X, X, X, X, X, X, X, X, X, X, X},{1, 3, 8, 9, 1, 8, X, X, X, X, X, X, X, X, X, X},{0, 9, 1, X, X, X, X, X, X, X, X, X, X, X, X, X},{0, 3, 8, X, X, X, X, X, X, X, X, X, X, X, X, X},{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}
};

边表用来获取在表面的顶点。而三角形表可以获得表面的三角片面序列。

Marcing Cube用以流体绘制

根据之前的SPH光滑核章节公式3.3:

我们计算得到的每个粒子的密度便可以作为标量场的值，该标量场也称为密度场。我们的粒子既然拥有自己的场，那么根据上面的方法便可以绘制流体的表面了。

MC类

class rxFace
{
public:vector<int> vert_idx;    // 顶点索引string material_name;    // 材质vector<glm::vec2> texcoords;    // 像素坐标int attribute;            //  属性public:rxFace() : attribute(0) {}inline int& operator[](int i){ return vert_idx[i]; }inline int  operator[](int i) const { return vert_idx[i]; }// 按函数进行顶点访问inline int& at(int i){ return vert_idx.at(i); }inline int  at(int i) const { return vert_idx.at(i); }// 更改多边形顶点的数量void resize(int size){vert_idx.resize(size);}// 返回多边形顶点的数量int size(void) const{return (int)vert_idx.size();}// 初始化void clear(void){vert_idx.clear();material_name = "";texcoords.clear();}
};struct RxVertexID{ //点的信息unsigned int newID;//索引double x,y,z;//位置
};typedef std::map<unsigned int,RxVertexID> ID2VertexID;struct RxTriangle{ //三角面的信息unsigned int pointID[3];//存有每个点的索引
};typedef std::vector<RxTriangle> RxTriangleVector;struct RxScalarField{//标量场unsigned int iNum[3];//x,y,z的大小glm::vec3 fWidth;//每个网格宽度glm::vec3 fMin;//最小点坐标
};class rxMCMesh{
public:rxMCMesh();~rxMCMesh();//从样本量生成三角形网格bool CreateMeshV(float *field, glm::vec3 min_p, double h, int n[3], float threshold,vector<glm::vec3> &vrts, vector<glm::vec3> &nrms, vector<rxFace> &face);//从样本量生成等值面网格void GenerateSurfaceV(const RxScalarField sf, float *field, float threshold,vector<glm::vec3> &vrts, vector<glm::vec3> &nrms, vector<int> &tris);//等值面创建成功则返回truebool IsSurfaceValid() const { return m_bValidSurface; }//删除表面void DeleteSurface();//用于网格划分的网格大小（在未创建网格的情况下，返回值为-1）int GetVolumeLengths(double& fVolLengthX, double& fVolLengthY, double& fVolLengthZ);//有关创建的网格的信息unsigned int GetNumVertices(void) const { return m_nVertices; }unsigned int GetNumTriangles(void) const { return m_nTriangles; }unsigned int GetNumNormals(void) const { return m_nNormals; }private:unsigned int m_nVertices;    //点的个数unsigned int m_nNormals;    //法线个数unsigned int m_nTriangles;    //三角面个数ID2VertexID m_i2pt3idVertices;            //形成等值面的顶点列表(以边索引为key，等值面的点为value)RxTriangleVector m_trivecTriangles;       //形成三角形的顶点列表RxScalarField m_Grid;            //标量场信息const float * m_ptScalarField;     //保存标量值的样本量float m_tIsoLevel;   //阈值bool m_bValidSurface;            //表面是否生成成功// 边idunsigned int GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo);// 顶点IDunsigned int GetVertexID(unsigned int nX, unsigned int nY, unsigned int nZ);// 计算边缘上的等值面RxVertexID CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nZ,unsigned int nEdgeNo);// 通过网格边缘两端的隐式函数值的线性插值计算等值点RxVertexID Interpolate(double fX1, double fY1, double fZ1, double fX2, double fY2, double fZ2,float tVal1, float tVal2);// 以输出形式存储顶点和网格几何信息void RenameVerticesAndTriangles(vector<glm::vec3> &vrts, unsigned int &nvrts, vector<int> &tris,unsigned int &ntris);// 顶点法线计算void CalculateNormals(const vector<glm::vec3> &vrts, unsigned int nvrts, const vector<int> &tris,unsigned int ntris,vector<glm::vec3> &nrms, unsigned int &nnrms);};

我们首先定义一个rxFace结构体，来保存每个面内的信息。定义RxVertexID作为保存点的信息。定义ID2VertexID为以网格边做为索引，表面点信息作为值的map。定义了RxTriangle的三角片面结构体，用以保存该三角片面中三个顶点索引。RxTriangleVector则用vector保存了所有的三角片面。rxMCMesh类中便有以ID2VertexID 定义的类成员，与RxTriangleVector定义的类成员。

由于我们绘制三角形的时候，是先根据边表，判断点落在哪条边上，然后通过插值计算出具体位置，然后我们可以通过ID2VertexID类对象保存以当前边为索引，索引的值为当前点点具体位置。而由于我们的三角形表获取的信息是点在哪条边上，所以我们用RxTriangleVector类对象从三角形表获取的边信息，加上通过ID2VertexID类对象从边信息获取的点具体位置，便可以获得最终三角形片面的三个顶点的最终位置了。

我们主要观察函数实现：

rxMCMesh::rxMCMesh(){m_Grid.fMin=glm::vec3(0.0);m_Grid.fWidth = glm::vec3(0.0);m_Grid.iNum[0] = 0;m_Grid.iNum[1] = 0;m_Grid.iNum[2] = 0;m_nTriangles = 0;m_nNormals = 0;m_nVertices = 0;m_ptScalarField = nullptr;m_fpScalarFunc = 0;m_tIsoLevel = 0;m_bValidSurface = false;
}rxMCMesh::~rxMCMesh(){DeleteSurface();
}void rxMCMesh::DeleteSurface()
{m_Grid.fWidth[0] = 0;m_Grid.fWidth[1] = 0;m_Grid.fWidth[2] = 0;m_Grid.iNum[0] = 0;m_Grid.iNum[1] = 0;m_Grid.iNum[2] = 0;m_nTriangles = 0;m_nNormals = 0;m_nVertices = 0;m_ptScalarField = NULL;m_tIsoLevel = 0;m_bValidSurface = false;
}//*从样本成三角形网格
//* @param [in]字段样本量
//* @param [in] min_p网格的最小坐标
//* @param [in] h网格的宽度
//* @param [in] n [3]网格编号（x，y，z）
//* @param [in]阈值阈值（隐含函数值网格化此值）
//* @param [in]方法网格生成方法（“mc”，“rmt”，“bloomenthal”）
//* @param [out] vrts顶点坐标
//* @param [out] nrms顶点正常
//* @param [out] tris mesh
//* @ retval true Mesh代成功
//* @retval false网格生成失败bool rxMCMesh::CreateMeshV(float *field, glm::vec3 min_p, double h, int *n, float threshold, vector<glm::vec3> &vrts, vector<glm::vec3> &nrms, vector<rxFace> &face){if(field == nullptr) return false;RxScalarField sf;for(int i = 0; i < 3; ++i){sf.iNum[i] = n[i];sf.fWidth[i] = h;sf.fMin[i] = min_p[i];}vector<int> tris;GenerateSurfaceV(sf, field, threshold, vrts, nrms, tris);if(IsSurfaceValid()){int nv = (int)GetNumVertices();int nm = (int)GetNumTriangles();int nn = (int)GetNumNormals();for(int i = 0; i < nn; ++i){nrms[i] *= -1.0;}face.resize(nm);for(int i = 0; i < nm; ++i){face[i].vert_idx.resize(3);for(int j = 0; j < 3; ++j){face[i][j] = tris[3*i+(2-j)];}}return true;}return false;
}void rxMCMesh::GenerateSurfaceV(const RxScalarField sf, float *field, float threshold,vector<glm::vec3> &vrts, vector<glm::vec3> &nrms, vector<int> &tris){// 等值面生成if(m_bValidSurface){DeleteSurface();}m_tIsoLevel = threshold;m_Grid.iNum[0] = sf.iNum[0];m_Grid.iNum[1] = sf.iNum[1];m_Grid.iNum[2] = sf.iNum[2];m_Grid.fWidth = sf.fWidth;m_Grid.fMin = sf.fMin;m_ptScalarField = field;unsigned int slice0 = (m_Grid.iNum[0] + 1);unsigned int slice1 = slice0*(m_Grid.iNum[1] + 1);// 等値面の生成for(unsigned int z = 0; z < m_Grid.iNum[2]; ++z){for(unsigned int y = 0; y < m_Grid.iNum[1]; ++y){for(unsigned int x = 0; x < m_Grid.iNum[0]; ++x){// 计算网格中的顶点放置信息表参考索引unsigned int tableIndex = 0;if(m_ptScalarField[z*slice1 + y*slice0 + x] < m_tIsoLevel)tableIndex |= 1;if(m_ptScalarField[z*slice1 + (y+1)*slice0 + x] < m_tIsoLevel)tableIndex |= 2;if(m_ptScalarField[z*slice1 + (y+1)*slice0 + (x+1)] < m_tIsoLevel)tableIndex |= 4;if(m_ptScalarField[z*slice1 + y*slice0 + (x+1)] < m_tIsoLevel)tableIndex |= 8;if(m_ptScalarField[(z+1)*slice1 + y*slice0 + x] < m_tIsoLevel)tableIndex |= 16;if(m_ptScalarField[(z+1)*slice1 + (y+1)*slice0 + x] < m_tIsoLevel)tableIndex |= 32;if(m_ptScalarField[(z+1)*slice1 + (y+1)*slice0 + (x+1)] < m_tIsoLevel)tableIndex |= 64;if(m_ptScalarField[(z+1)*slice1 + y*slice0 + (x+1)] < m_tIsoLevel)tableIndex |= 128;if(edgeTable[tableIndex] != 0){// 计算边上的顶点if(edgeTable[tableIndex] & 8){RxVertexID pt = CalculateIntersection(x, y, z, 3);unsigned int id = GetEdgeID(x, y, z, 3);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if(edgeTable[tableIndex] & 1){RxVertexID pt = CalculateIntersection(x, y, z, 0);unsigned int id = GetEdgeID(x, y, z, 0);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if(edgeTable[tableIndex] & 256){RxVertexID pt = CalculateIntersection(x, y, z, 8);unsigned int id = GetEdgeID(x, y, z, 8);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if(x == m_Grid.iNum[0] - 1){if(edgeTable[tableIndex] & 4){RxVertexID pt = CalculateIntersection(x, y, z, 2);unsigned int id = GetEdgeID(x, y, z, 2);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if(edgeTable[tableIndex] & 2048){RxVertexID pt = CalculateIntersection(x, y, z, 11);unsigned int id = GetEdgeID(x, y, z, 11);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}}if(y == m_Grid.iNum[1] - 1){if(edgeTable[tableIndex] & 2){RxVertexID pt = CalculateIntersection(x, y, z, 1);unsigned int id = GetEdgeID(x, y, z, 1);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if(edgeTable[tableIndex] & 512){RxVertexID pt = CalculateIntersection(x, y, z, 9);unsigned int id = GetEdgeID(x, y, z, 9);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}}if(z == m_Grid.iNum[2] - 1){if(edgeTable[tableIndex] & 16){RxVertexID pt = CalculateIntersection(x, y, z, 4);unsigned int id = GetEdgeID(x, y, z, 4);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if(edgeTable[tableIndex] & 128){RxVertexID pt = CalculateIntersection(x, y, z, 7);unsigned int id = GetEdgeID(x, y, z, 7);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}}if((x==m_Grid.iNum[0] - 1) && (y==m_Grid.iNum[1] - 1))if(edgeTable[tableIndex] & 1024){RxVertexID pt = CalculateIntersection(x, y, z, 10);unsigned int id = GetEdgeID(x, y, z, 10);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if((x==m_Grid.iNum[0] - 1) && (z==m_Grid.iNum[2] - 1))if(edgeTable[tableIndex] & 64){RxVertexID pt = CalculateIntersection(x, y, z, 6);unsigned int id = GetEdgeID(x, y, z, 6);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}if((y==m_Grid.iNum[1] - 1) && (z==m_Grid.iNum[2] - 1))if(edgeTable[tableIndex] & 32){RxVertexID pt = CalculateIntersection(x, y, z, 5);unsigned int id = GetEdgeID(x, y, z, 5);m_i2pt3idVertices.insert(ID2VertexID::value_type(id, pt));}// 多边形生成for(unsigned int i = 0; triTable[tableIndex][i] != 255; i += 3){RxTriangle triangle;unsigned int pointID0, pointID1, pointID2;pointID0 = GetEdgeID(x, y, z, triTable[tableIndex][i]);pointID1 = GetEdgeID(x, y, z, triTable[tableIndex][i+1]);pointID2 = GetEdgeID(x, y, z, triTable[tableIndex][i+2]);triangle.pointID[0] = pointID0;triangle.pointID[1] = pointID1;triangle.pointID[2] = pointID2;m_trivecTriangles.push_back(triangle);}}}}}RenameVerticesAndTriangles(vrts, m_nVertices, tris, m_nTriangles);CalculateNormals(vrts, m_nVertices, tris, m_nTriangles, nrms, m_nNormals);m_bValidSurface = true;
}/**获得边缘ID* @param [in] nX，nY，nZ网格位置* @param [in] nEdgeNo边数* @return边缘ID*/
unsigned int rxMCMesh::GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
{switch(nEdgeNo){case 0:return GetVertexID(nX, nY, nZ) + 1;case 1:return GetVertexID(nX, nY + 1, nZ);case 2:return GetVertexID(nX + 1, nY, nZ) + 1;case 3:return GetVertexID(nX, nY, nZ);case 4:return GetVertexID(nX, nY, nZ + 1) + 1;case 5:return GetVertexID(nX, nY + 1, nZ + 1);case 6:return GetVertexID(nX + 1, nY, nZ + 1) + 1;case 7:return GetVertexID(nX, nY, nZ + 1);case 8:return GetVertexID(nX, nY, nZ) + 2;case 9:return GetVertexID(nX, nY + 1, nZ) + 2;case 10:return GetVertexID(nX + 1, nY + 1, nZ) + 2;case 11:return GetVertexID(nX + 1, nY, nZ) + 2;default:// Invalid edge no.return -1;}
}/*！
*获取顶点ID
* @param [in] nX，nY，nZ网格位置
* @return顶点ID
*/
unsigned int rxMCMesh::GetVertexID(unsigned int nX, unsigned int nY, unsigned int nZ)
{return 3*(nZ*(m_Grid.iNum[1] + 1)*(m_Grid.iNum[0] + 1) + nY*(m_Grid.iNum[0] + 1) + nX);
}/*！* 通过插值计算边缘上的等值点（从样本量）* @param [in] nX，nY，nZ网格位置* @param [in] nEdgeNo边数* @return网格顶点信息*/
RxVertexID rxMCMesh::CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
{double x1, y1, z1, x2, y2, z2;unsigned int v1x = nX, v1y = nY, v1z = nZ;unsigned int v2x = nX, v2y = nY, v2z = nZ;switch(nEdgeNo){case 0:v2y += 1;break;case 1:v1y += 1;v2x += 1;v2y += 1;break;case 2:v1x += 1;v1y += 1;v2x += 1;break;case 3:v1x += 1;break;case 4:v1z += 1;v2y += 1;v2z += 1;break;case 5:v1y += 1;v1z += 1;v2x += 1;v2y += 1;v2z += 1;break;case 6:v1x += 1;v1y += 1;v1z += 1;v2x += 1;v2z += 1;break;case 7:v1x += 1;v1z += 1;v2z += 1;break;case 8:v2z += 1;break;case 9:v1y += 1;v2y += 1;v2z += 1;break;case 10:v1x += 1;v1y += 1;v2x += 1;v2y += 1;v2z += 1;break;case 11:v1x += 1;v2x += 1;v2z += 1;break;}// 获取边的两点坐标x1 = m_Grid.fMin[0]+v1x*m_Grid.fWidth[0];y1 = m_Grid.fMin[1]+v1y*m_Grid.fWidth[1];z1 = m_Grid.fMin[2]+v1z*m_Grid.fWidth[2];x2 = m_Grid.fMin[0]+v2x*m_Grid.fWidth[0];y2 = m_Grid.fMin[1]+v2y*m_Grid.fWidth[1];z2 = m_Grid.fMin[2]+v2z*m_Grid.fWidth[2];unsigned int slice0 = (m_Grid.iNum[0] + 1);unsigned int slice1 = slice0*(m_Grid.iNum[1] + 1);float val1 = m_ptScalarField[v1z*slice1 + v1y*slice0 + v1x];float val2 = m_ptScalarField[v2z*slice1 + v2y*slice0 + v2x];RxVertexID intersection = Interpolate(x1, y1, z1, x2, y2, z2, val1, val2);return intersection;
}/*！*网格顶点，以输出格式存储的几何信息* @param [out] vrts顶点坐标* @param [out] nvrts顶点数* @param [out] tris三角形多边形几何信息* @param [out] ntris三角形多边形的数量*/
void rxMCMesh::RenameVerticesAndTriangles(vector<glm::vec3> &vrts, unsigned int &nvrts, vector<int> &tris, unsigned int &ntris)
{unsigned int nextID = 0;ID2VertexID::iterator mapIterator = m_i2pt3idVertices.begin();RxTriangleVector::iterator vecIterator = m_trivecTriangles.begin();// 刷新点while(mapIterator != m_i2pt3idVertices.end()){(*mapIterator).second.newID = nextID;nextID++;mapIterator++;}// 刷新三角面.while(vecIterator != m_trivecTriangles.end()){for(unsigned int i = 0; i < 3; i++){unsigned int newID = m_i2pt3idVertices[(*vecIterator).pointID[i]].newID;(*vecIterator).pointID[i] = newID;}vecIterator++;}// 将所有顶点和三角形复制到两个数组中，以便可以有效地访问它们。// 复制点mapIterator = m_i2pt3idVertices.begin();nvrts = (int)m_i2pt3idVertices.size();vrts.resize(nvrts);for(unsigned int i = 0; i < nvrts; i++, mapIterator++){vrts[i][0] = (*mapIterator).second.x;vrts[i][1] = (*mapIterator).second.y;vrts[i][2] = (*mapIterator).second.z;}// 复制制作三角形的顶点索引。vecIterator = m_trivecTriangles.begin();ntris = (int)m_trivecTriangles.size();tris.resize(ntris*3);for(unsigned int i = 0; i < ntris; i++, vecIterator++){tris[3*i+0] = (*vecIterator).pointID[0];tris[3*i+1] = (*vecIterator).pointID[1];tris[3*i+2] = (*vecIterator).pointID[2];}//释放空间m_i2pt3idVertices.clear();m_trivecTriangles.clear();
}/*！*顶点法线计算* @param [in] vrts顶点坐标* @param [in] nvrts顶点数* @param [in] tris三角形多边形几何信息* @param [in] ntris三角形多边形的数量* @param [out] nrms法线* @param [out] nnrms法线数（=顶点数）*/
void rxMCMesh::CalculateNormals(const vector<glm::vec3> &vrts, unsigned int nvrts, const vector<int> &tris, unsigned int ntris,vector<glm::vec3> &nrms, unsigned int &nnrms)
{nnrms = nvrts;nrms.resize(nnrms);// 初始化法线.for( unsigned int i = 0; i < nnrms; i++){nrms[i][0] = 0;nrms[i][1] = 0;nrms[i][2] = 0;}// 计算法线.for( unsigned int i = 0; i < ntris; i++){glm::vec3 vec1, vec2, normal;unsigned int id0, id1, id2;id0 = tris[3*i+0];id1 = tris[3*i+1];id2 = tris[3*i+2];vec1 = vrts[id1]-vrts[id0];vec2 = vrts[id2]-vrts[id0];normal = cross(vec1, vec2);nrms[id0] += normal;nrms[id1] += normal;nrms[id2] += normal;}// 单位化法线.for( unsigned int i = 0; i < nnrms; i++){normalize(nrms[i]);}
}

其中GenerateSurfaceV函数则是我们的Marching Cube的应用过程。我们把密度场作为输入，遍历密度场的每一个网格点，判断其是否超过密度阈值（超过标1否则标0）。最后的tableIndex的八位2进制则保存了顶点索引。然后我们将tableIndex作为边表edgeTable的索引，用它获取其在网格的哪一条边上并用CalculateIntersection函数计算出线性插值，确定其位置，并得到该点的对应边编号。我们最后将该边编号作为索引，点位置作为值传入m_i2pt3idVertices。另外我们还用三角形表索引出每个网格内的三角形序列，并保存至m_trivecTriangles。有了每个三角片面的每个顶点，我们就能对他们求取法线了，函数CalculateNormals实现了该过程。

之后我们如果只要表面的顶点，只需要访问m_i2pt3idVertices内的所有点。若需要构建表面，只需要访问m_trivecTriangles即可。

到这里我们的MC类就全完成了，接下来的一小节我们只要更改一下我们之前的流体系统，为其添加mc的应用，那么我们的流体就不会再是一个个点，而是有自己的表面了。

SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（1）相关推荐

SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（2）
流体模拟(三) Marching Cube算法(2) 我们在之前的流体系统类里新加入一些函数和成员,用来引用我们的MC类,便可以获得生成有表面的流体模型了,效果如图: 添加简单的天空盒以及Phong氏 ...
SPH（光滑粒子流体动力学）流体模拟实现六：Position Based Fluid（PBF）
SPH(光滑粒子流体动力学)流体模拟实现六:Position Based Fluid(PBF) PBF方法和前篇提到的PCISPH方法类似,都属于迭代矫正法.PCISPH是通过迭代预测压力,通过压力变 ...
SPH（光滑粒子流体动力学）流体模拟实现七：屏幕空间流体渲染（SSF）
SPH(光滑粒子流体动力学)流体模拟实现七:屏幕空间流体渲染(SSF) 之前都是用的Marching Cube重建流体表面的方法.近来为了做对比实验,尝试了屏幕空间流体渲染的方法.发现屏幕空间的方法不 ...
SPH（光滑粒子流体动力学）流体模拟实现：算法总览
流体模拟(一) 流体模拟算法总体流程: 流体现象广泛存在于自然界.日常生活以及工业生产中,对流体的模拟即流体动画, 一直是基于物理的动画以及计算机图形学的重要研究内容.目前, 基于物理模拟的流体动画 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现2
流体模拟(二) SPH算法实现2: 在前面一节我们完成了粒子缓存类,网格类和邻接表类.我们现在可以正式的整合在我们的流体系统类中了. 流体系统类 class FluidSystem{public:Fl ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现1
流体模拟(二) SPH算法实现1: 由于我们计算每个粒子的状态时,都需要获得在它光滑核半径内(邻域内)的所有粒子信息.我们如果遍历每个粒子,计算欧式距离的话,那开销就过于庞大了.因此我们可以将我们的空 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（3）-光滑核函数
流体模拟(二) 光滑核函数: sph中涉及的光滑核可以理解为:在一定的光滑核半径内,所受的力受距离权重的影响,距离越近所受影响越大.其表现形式如图所示. 这里我们便可以将流体看成一个个粒子的集合,每一 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（2）-粒子受力分析
流体模拟(二) SPH算法的粒子受力分析: SPH算法的基本设想,就是将连续的流体想象成一个个相互作用的微粒,这些粒子相互影响,共同形成了复杂的流体运动.其实现的原理则是我们在初始空间里创建多个粒子, ...
SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(1)
流体模拟四: 各向异性(Anisotropic)表面光滑(1) 表面的表示与定义我们是以隐式表面来表示流体表面.定义一个足够大的参考域,它完全包含粒子表示的流体所在的空间D,即.流体表面隐含的定义为 ...

SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（1）

流体模拟（三）

Marching Cube算法（1）

Marcing Cube用以流体绘制

MC类

SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（1）相关推荐

最新文章

热门文章