流体模拟(三)
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个顶点,那所有的可能性就只有=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);
//等值面创建成功则返回true
bool 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; //表面是否生成成功
// 边id
unsigned int GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo);
// 顶点ID
unsigned 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的应用,那么我们的流体就不会再是一个个点,而是有自己的表面了。