Marching Cubes 算法再探
- 基础
- 过程
- 代码
- mian.cpp
- MarchingCubes.h
- MarchingCubes.cpp
之前做NeRF相关工作时简单看过，但没有深究其实现，知其然不知其所以然的程度，算是初探。
基础
为了对MC有大致了解，可以先根据Marching Cubes 算法初探，理解一下Marching Cubes这个名字的由来，其二维空间中的示例可谓一目了然。
接下来结合Efficient Implementation of Marching Cubes’ Cases with Topological Guarantees
的论文和代码，来探究一下CPU版本经典MC的具体实现过程。
过程
整个过程的核心部分可以抽象为以下步骤：
- 对每一个cube，比较其顶点与等值面的大小，得到顶点的8bit案例编码
顶点1、5、7值为正，则得到10100010，即162
- 根据顶点案例的8bit编码在查找表中找到相应三角形数组
162对应的三角形数组为{5,0,1,5,4,0,7,6,11}
- 根据三角形数组添加\绘制三角形
这个三角形数组每三个值代表一个三角形的三个顶点位置（先以cube边索引形式给出）
添加三角形[5,0,1]
添加三角形[5,4,0]
添加三角形[7,6,11]
- 遍历所有cube，得到所有三角形的并集
以上图截取于视频，这个视频也很有意思。
代码
接下来看下大佬Thomas Lewiner(Efficient Implementation of Marching Cubes’ Cases with Topological Guarantees一作)的开源代码。
mian.cpp
//_____________________________________________________________________________
// main function
int main (int argc, char **argv)
//-----------------------------------------------------------------------------
{
MarchingCubes mc ; // 构造
mc.set_resolution( 60,60,60 ) ; // 设置分辨率
//mc.set_method(true); // 是否使用原始MC方法
mc.init_all() ; // 初始化
compute_data( mc ) ; // 生成模拟数据，调用 mc.set_data()
mc.run() ; // 主函数
mc.clean_temps() ; // 清除临时变量
mc.writePLY("test.ply") ; // 保存结果
mc.clean_all() ; // 析构
return 0 ;
}
//_____________________________________________________________________________
这位大佬写的代码相当规范且有非常详尽的注释，考虑篇幅，在下面仅展示关键部分
MarchingCubes.h
typedef unsigned char uchar ;
typedef signed char schar ;
typedef float real ;
class MarchingCubes
//-----------------------------------------------------------------------------
{
// 构造函数
public :
/** 主构造函数和默认构造函数*/
MarchingCubes ( const int size_x = -1, const int size_y = -1, const int size_z = -1 ) ;
/** 析构函数 */
~MarchingCubes() ;
//-----------------------------------------------------------------------------
// 访问器
public :
/**
* 修改网格的大小
* \param size_x 网格的宽度
* \param size_y 网格的深度
* \param size_z 网格的高度
*/
inline void set_resolution( const int size_x, const int size_y, const int size_z )
{
_size_x = size_x ;
_size_y = size_y ;
_size_z = size_z ;
}
/**
* 选择是否使用增强的拓扑控制查找表或原始的 Marching Cubes 算法
* \param originalMC 如果为 true，则使用原始的 Marching Cubes
*/
inline void set_method ( const bool originalMC = false ) { _originalMC = originalMC ; }
/**
* 设置网格中指定立方体的数据
* \param val 新的立方体值
* \param i 立方体的横坐标
* \param j 立方体的纵坐标
* \param k 立方体的高度
*/
inline void set_data ( const real val, const int i, const int j, const int k )
{
_data[ i + j*_size_x + k*_size_x*_size_y] = val ;
}
// 数据初始化
/** 初始化临时结构（必须在调用前设置大小）：包括网格和每个立方体的顶点索引 */
void init_temps();
/** 初始化所有结构（必须在调用前设置大小）：包括临时结构和网格缓存 */
void init_all();
/** 清除临时结构：包括网格和主要结构 */
void clean_temps();
/** 清除所有结构：包括临时结构和网格缓存 */
void clean_all();
//-----------------------------------------------------------------------------
// 算法
public :
/**
* 主算法：必须在调用 init_all 之后调用
* \param iso 等值面值
*/
void run( real iso = (real)0.0 ) ;
protected :
/** 处理一个立方体的三角剖分 */
void process_cube () ;
/** 测试立方体的三角剖分组件是否应通过模糊面的内部连接 */
bool test_face ( schar face ) ;
/** 测试立方体的三角剖分组件是否应通过立方体的内部连接 */
bool test_interior( schar s ) ;
//-----------------------------------------------------------------------------
// 操作
protected :
/**
* 通过沿立方体边缘的插值计算几乎所有网格的顶点
* \param iso 等值面值
*/
void compute_intersection_points( real iso ) ;
/**
* 添加三角形到网格的例程
* \param trig 三角形的代码作为边索引的序列
* \param n 需要生成的三角形数量
* \param v12 使用的内部顶点的索引（如有必要）
*/
void add_triangle ( const char* trig, char n, int v12 = -1 ) ;
/** 测试并在需要时增加顶点缓存的容量以插入新顶点 */
void test_vertex_addition() ;
/** 在当前水平边上添加顶点 */
int add_x_vertex() ;
//yz
/** 在当前立方体内添加顶点 */
int add_c_vertex() ;
/**
* 插值指定立方体下顶点的隐函数的水平梯度
* \param i 立方体的横坐标
* \param j 立方体的纵坐标
* \param k 立方体的高度
*/
real get_x_grad( const int i, const int j, const int k ) const ;
//yz
/**
* 获取指定立方体下水平边的预计算顶点索引
* \param i 立方体的横坐标
* \param j 立方体的纵坐标
* \param k 立方体的高度
*/
inline int get_x_vert( const int i, const int j, const int k ) const
{
return _x_verts[ i + j*_size_x + k*_size_x*_size_y] ;
}
//yz
/**
* 设置特定立方体下水平边的预计算顶点索引
* \param val 新顶点的索引
* \param i 立方体的横坐标
* \param j 立方体的纵坐标（平面内的纵向）
* \param k 立方体的高度（深度方向的坐标）
*/
inline void set_x_vert(const int val, const int i, const int j, const int k)
{
_x_verts[i + j * _size_x + k * _size_x * _size_y] = val;
}
//yz
//-----------------------------------------------------------------------------
// 元素
protected :
bool _originalMC ; /**< 选择算法是使用增强的拓扑控制查找表还是使用原始的 Marching Cubes 算法 */
bool _ext_data ; /**< 选择是分配数据还是使用来自其他类的数据 */
int _size_x ; /**< 网格的宽度 */
int _size_y ; /**< 网格的深度 */
int _size_z ; /**< 网格的高度 */
real *_data ; /**< 在网格上采样的隐函数值 */
int *_x_verts ; /**< 预计算的每个立方体下水平边的顶点索引 */
int *_y_verts ; /**< 预计算的每个立方体下纵向边的顶点索引 */
int *_z_verts ; /**< 预计算的每个立方体下垂直边的顶点索引 */
int _nverts ; /**< 顶点缓存中已分配的顶点数量 */
int _ntrigs ; /**< 三角形缓存中已分配的三角形数量 */
int _Nverts ; /**< 顶点缓存的大小 */
int _Ntrigs ; /**< 三角形缓存的大小 */
Vertex *_vertices ; /**< 顶点缓存 */
Triangle *_triangles ; /**< 三角形缓存 */
int _i ; /**< 当前活动立方体的横坐标-对应x */
int _j ; /**< 当前活动立方体的高度-对应y */
int _k ; /**< 当前活动立方体的纵坐标-对应z */
real _cube[8] ; /**< 当前活动立方体上的隐函数值 */
uchar _lut_entry ; /**< 立方体的符号表示，范围在 [0..255] */
uchar _case ; /**< 当前活动立方体的案例，范围在 [0..15] */
uchar _config ; /**< 当前活动立方体的配置 */
uchar _subconfig ; /**< 当前活动立方体的子配置 */
};
MarchingCubes.cpp
//-----------------------------------------------------------------------------
// 初始化临时结构（调用前需要设置大小）
void MarchingCubes::init_temps()
//-----------------------------------------------------------------------------
{
if( !_ext_data ) // 如果没有使用外部数据
_data = new real [_size_x * _size_y * _size_z] ; // 分配网格数据的内存
_x_verts = new int [_size_x * _size_y * _size_z] ; // 分配用于存储水平边顶点索引的数组
_y_verts = new int [_size_x * _size_y * _size_z] ; // 分配用于存储纵向边顶点索引的数组
_z_verts = new int [_size_x * _size_y * _size_z] ; // 分配用于存储垂直边顶点索引的数组
memset( _x_verts, -1, _size_x * _size_y * _size_z * sizeof( int ) ) ; // 初始化水平边顶点索引数组
memset( _y_verts, -1, _size_x * _size_y * _size_z * sizeof( int ) ) ; // 初始化纵向边顶点索引数组
memset( _z_verts, -1, _size_x * _size_y * _size_z * sizeof( int ) ) ; // 初始化垂直边顶点索引数组
}
//_____________________________________________________________________________
//_____________________________________________________________________________
// 初始化所有结构（调用前需要设置大小）
void MarchingCubes::init_all ()
//-----------------------------------------------------------------------------
{
init_temps() ; // 初始化临时结构
_nverts = _ntrigs = 0 ; // 初始化顶点数和三角形数为 0
_Nverts = _Ntrigs = ALLOC_SIZE ; // 设置顶点缓冲区和三角形缓冲区的初始大小(#define ALLOC_SIZE 65536)
_vertices = new Vertex [_Nverts] ; // 分配顶点缓冲区的内存
_triangles = new Triangle[_Ntrigs] ; // 分配三角形缓冲区的内存
}
//_____________________________________________________________________________
在以下函数中，根据等效旧代码可以明显看出for循环内是在确定每个cube的8bit编码(已经转为十进制)即所谓的查找表入口(_lut_entry)，再根据编码找到相应的案例，再经过拓扑验证等操作得到每个cube的三角近似
//_____________________________________________________________________________
// 主算法
void MarchingCubes::run( real iso )
//-----------------------------------------------------------------------------
{
clock_t time = clock() ; // 记录开始时间
compute_intersection_points( iso ) ; // 计算网格中每个立方体边界上的交点
// 遍历三维网格中的每个立方体
for( _k = 0 ; _k < _size_z-1 ; _k++ )
for( _j = 0 ; _j < _size_y-1 ; _j++ )
for( _i = 0 ; _i < _size_x-1 ; _i++ )
{
_lut_entry = 0 ; // 初始化查找表的入口值
for( int p = 0 ; p < 8 ; ++p )
{
// 计算当前立方体的每个顶点相对于等值面的值
_cube[p] = get_data( _i+((p^(p>>1))&1), _j+((p>>1)&1), _k+((p>>2)&1) ) - iso ;
if( fabs( _cube[p] ) < FLT_EPSILON ) _cube[p] = FLT_EPSILON ; // 如果值接近零，设置为最小浮点值
if( _cube[p] > 0 ) _lut_entry += 1 << p ; // 根据顶点值更新查找表入口
}
/*
// 以下是等效的旧代码，每个顶点的计算逐一展开
if( ( _cube[0] = get_data( _i , _j , _k ) ) > 0 ) _lut_entry += 1 ;
if( ( _cube[1] = get_data(_i+1, _j , _k ) ) > 0 ) _lut_entry += 2 ;
if( ( _cube[2] = get_data(_i+1,_j+1, _k ) ) > 0 ) _lut_entry += 4 ;
if( ( _cube[3] = get_data( _i ,_j+1, _k ) ) > 0 ) _lut_entry += 8 ;
if( ( _cube[4] = get_data( _i , _j ,_k+1) ) > 0 ) _lut_entry += 16 ;
if( ( _cube[5] = get_data(_i+1, _j ,_k+1) ) > 0 ) _lut_entry += 32 ;
if( ( _cube[6] = get_data(_i+1,_j+1,_k+1) ) > 0 ) _lut_entry += 64 ;
if( ( _cube[7] = get_data( _i ,_j+1,_k+1) ) > 0 ) _lut_entry += 128 ;
*/
process_cube() ; // 处理当前立方体，生成对应的三角面片
}
// 打印算法运行时间
printf("Marching Cubes ran in %lf secs.\n", (double)(clock() - time)/CLOCKS_PER_SEC) ;
}
//_____________________________________________________________________________
根据上面顶点的编码(第一张截图)，就可以理解为什么以下过程中只处理顶点0、1、3、4
下面的代码的大意就是，循环每个处理cube的单位轴，但凡轴边两端点正负值不同，则在相应边上添加顶点
//-----------------------------------------------------------------------------
// 计算等值面交点
void MarchingCubes::compute_intersection_points( real iso )
//-----------------------------------------------------------------------------
{
for( _k = 0 ; _k < _size_z ; _k++ ) // 遍历 z 方向
for( _j = 0 ; _j < _size_y ; _j++ ) // 遍历 y 方向
for( _i = 0 ; _i < _size_x ; _i++ ) // 遍历 x 方向
{
_cube[0] = get_data( _i, _j, _k ) - iso ; // 获取当前顶点的值并减去等值面值
if( _i < _size_x - 1 ) // 如果不是x方向的最后一个cube
_cube[1] = get_data(_i+1, _j , _k ) - iso ; // 获取下一顶点的值并减去等值面值
else // x方向的最后一个cube
_cube[1] = _cube[0] ; // 使用上一顶点的值
if( _j < _size_y - 1 ) // 如果不是y方向的最后一个cube
_cube[3] = get_data( _i ,_j+1, _k ) - iso ; // 获取下一顶点的值并减去等值面值
else // 如果是
_cube[3] = _cube[0] ; // 使用上一顶点的值
if( _k < _size_z - 1 ) // 如果不是z方向的最后一个cube
_cube[4] = get_data( _i , _j ,_k+1) - iso ; // 获取下一顶点的值并减去等值面值
else // 如果是
_cube[4] = _cube[0] ; // 使用上一顶点的值
if( fabs( _cube[0] ) < FLT_EPSILON ) _cube[0] = FLT_EPSILON ; // 如果接近于零，则设置为最小正数
if( fabs( _cube[1] ) < FLT_EPSILON ) _cube[1] = FLT_EPSILON ;
if( fabs( _cube[3] ) < FLT_EPSILON ) _cube[3] = FLT_EPSILON ;
if( fabs( _cube[4] ) < FLT_EPSILON ) _cube[4] = FLT_EPSILON ;
if( _cube[0] < 0 ) // 如果当前顶点的值小于等值面值
{
// 如果三个方向的下一顶点的值大于等值面值，则在相应边上添加顶点
if( _cube[1] > 0 ) set_x_vert( add_x_vertex( ), _i,_j,_k ) ;
if( _cube[3] > 0 ) set_y_vert( add_y_vertex( ), _i,_j,_k ) ;
if( _cube[4] > 0 ) set_z_vert( add_z_vertex( ), _i,_j,_k ) ;
}
else // 如果当前顶点的值大于等于等值面值
{
// 如果三个方向的下一顶点的值小于等值面值，则在相应边上添加顶点
if( _cube[1] < 0 ) set_x_vert( add_x_vertex( ), _i,_j,_k ) ;
if( _cube[3] < 0 ) set_y_vert( add_y_vertex( ), _i,_j,_k ) ;
if( _cube[4] < 0 ) set_z_vert( add_z_vertex( ), _i,_j,_k ) ;
}
}
}
//_____________________________________________________________________________
根据案例编码确定案例，再通过配置（复杂的案例需要子配置）保证正确的拓扑，部分代码太长了被省略，这里也不一一看各种案例的具体情况了
//_____________________________________________________________________________
// 处理单位立方体
void MarchingCubes::process_cube( )
//-----------------------------------------------------------------------------
{
// 如果启用了经典的Marching Cubes模式
if( _originalMC )
{
char nt = 0 ;
// 计算当前查找表中三角形的数量
while( casesClassic[_lut_entry][3*nt] != -1 ) nt++ ;
// 添加三角形
add_triangle( casesClassic[_lut_entry], nt ) ;
return ;
}
int v12 = -1 ; // 用于存储生成的额外顶点索引
_case = cases[_lut_entry][0] ; // 获取当前立方体的拓扑情况（最大值14）
_config = cases[_lut_entry][1] ; // 获取当前立方体的配置（最大值47）
_subconfig = 0 ; // 初始化子配置
switch( _case )
{
case 0 : // 情况0：无需处理
break ;
case 1 : // 情况1：添加一个三角形
add_triangle( tiling1[_config], 1 ) ;
break ;
case 2 : // 情况2：添加两个三角形
add_triangle( tiling2[_config], 2 ) ;
break ;
case 3 : // 情况3：根据面测试的结果选择不同的三角形配置
if( test_face( test3[_config]) )
add_triangle( tiling3_2[_config], 4 ) ; // 3.2 配置
else
add_triangle( tiling3_1[_config], 2 ) ; // 3.1 配置
break ;
case 4 : // 情况4：根据内部测试的结果选择不同的三角形配置
if( test_interior( test4[_config]) )
add_triangle( tiling4_1[_config], 2 ) ; // 4.1.1 配置
else
add_triangle( tiling4_2[_config], 6 ) ; // 4.1.2 配置
break ;
case 5 :
add_triangle( tiling5[_config], 3 ) ;
break ;
case 6 :
if( test_face( test6[_config][0]) )
add_triangle( tiling6_2[_config], 5 ) ; // 6.2
else
{
if( test_interior( test6[_config][1]) )
add_triangle( tiling6_1_1[_config], 3 ) ; // 6.1.1
else
{
v12 = add_c_vertex() ;
add_triangle( tiling6_1_2[_config], 9 , v12) ; // 6.1.2
}
}
break ;
case 7 :
if( test_face( test7[_config][0] ) ) _subconfig += 1 ;
if( test_face( test7[_config][1] ) ) _subconfig += 2 ;
if( test_face( test7[_config][2] ) ) _subconfig += 4 ;
switch( _subconfig )
{
case 0 :
add_triangle( tiling7_1[_config], 3 ) ; break ;
// ......
case 7 :
if( test_interior( test7[_config][3]) )
add_triangle( tiling7_4_2[_config], 9 ) ;
else
add_triangle( tiling7_4_1[_config], 5 ) ;
break ;
};
break ;
case 8 :
add_triangle( tiling8[_config], 2 ) ;
break ;
case 9 :
add_triangle( tiling9[_config], 4 ) ;
break ;
case 10 :
if( test_face( test10[_config][0]) )
{
if( test_face( test10[_config][1]) )
add_triangle( tiling10_1_1_[_config], 4 ) ; // 10.1.1
else
{
v12 = add_c_vertex() ;
add_triangle( tiling10_2[_config], 8, v12 ) ; // 10.2
}
}
else
{
if( test_face( test10[_config][1]) )
{
v12 = add_c_vertex() ;
add_triangle( tiling10_2_[_config], 8, v12 ) ; // 10.2
}
else
{
if( test_interior( test10[_config][2]) )
add_triangle( tiling10_1_1[_config], 4 ) ; // 10.1.1
else
add_triangle( tiling10_1_2[_config], 8 ) ; // 10.1.2
}
}
break ;
case 11 :
add_triangle( tiling11[_config], 4 ) ;
break ;
case 12 :
if( test_face( test12[_config][0]) )
{
if( test_face( test12[_config][1]) )
add_triangle( tiling12_1_1_[_config], 4 ) ; // 12.1.1
else
{
v12 = add_c_vertex() ;
add_triangle( tiling12_2[_config], 8, v12 ) ; // 12.2
}
}
else
{
if( test_face( test12[_config][1]) )
{
v12 = add_c_vertex() ;
add_triangle( tiling12_2_[_config], 8, v12 ) ; // 12.2
}
else
{
if( test_interior( test12[_config][2]) )
add_triangle( tiling12_1_1[_config], 4 ) ; // 12.1.1
else
add_triangle( tiling12_1_2[_config], 8 ) ; // 12.1.2
}
}
break ;
case 13 :
if( test_face( test13[_config][0] ) ) _subconfig += 1 ;
if( test_face( test13[_config][1] ) ) _subconfig += 2 ;
if( test_face( test13[_config][2] ) ) _subconfig += 4 ;
if( test_face( test13[_config][3] ) ) _subconfig += 8 ;
if( test_face( test13[_config][4] ) ) _subconfig += 16 ;
if( test_face( test13[_config][5] ) ) _subconfig += 32 ;
switch( subconfig13[_subconfig] )
{
case 0 :/* 13.1 */
add_triangle( tiling13_1[_config], 4 ) ; break ;
// ......
case 45 :/* 13.1 */
add_triangle( tiling13_1_[_config], 4 ) ; break ;
default :
printf("Marching Cubes: Impossible case 13?\n" ) ; print_cube() ;
}
break ;
case 14 :
add_triangle( tiling14[_config], 4 ) ;
break ;
};
}
//_____________________________________________________________________________
test_face和test_interior都是用以确定正确的拓扑，具体思路可见原文
//-----------------------------------------------------------------------------
// 测试一个面
// 如果面大于0，则返回true 表示该面包含等值面的一部分
bool MarchingCubes::test_face( schar face )
//-----------------------------------------------------------------------------
{
real A, B, C, D;
switch( face )
{
case -1 : case 1 : A = _cube[0] ; B = _cube[4] ; C = _cube[5] ; D = _cube[1] ; break ; // x轴方向的面
case -2 : case 2 : A = _cube[1] ; B = _cube[5] ; C = _cube[6] ; D = _cube[2] ; break ; // y轴方向的面
case -3 : case 3 : A = _cube[2] ; B = _cube[6] ; C = _cube[7] ; D = _cube[3] ; break ; // z轴方向的面
case -4 : case 4 : A = _cube[3] ; B = _cube[7] ; C = _cube[4] ; D = _cube[0] ; break ; // 反x轴方向的面
case -5 : case 5 : A = _cube[0] ; B = _cube[3] ; C = _cube[2] ; D = _cube[1] ; break ; // 反y轴方向的面
case -6 : case 6 : A = _cube[4] ; B = _cube[7] ; C = _cube[6] ; D = _cube[5] ; break ; // 反z轴方向的面
default : printf( "Invalid face code %d\n", face ) ; print_cube() ; A = B = C = D = 0 ;
};
if( fabs( A*C - B*D ) < FLT_EPSILON ) // 如果AC-BD接近于零
return face >= 0 ; // 直接返回face的正负号
return face * A * ( A*C - B*D ) >= 0 ; // 如果face和A的符号相同，并且AC-BD大于零，则返回true
}
//_____________________________________________________________________________
//_____________________________________________________________________________
// 测试立方体的内部
// 如果 s == 7，则在内部为空时返回 true
// 如果 s == -7，则在内部为空时返回 false
bool MarchingCubes::test_interior(schar s)
//-----------------------------------------------------------------------------
{
real t, At=0, Bt=0, Ct=0, Dt=0, a, b;
char test = 0;
char edge = -1; // 三角剖分的参考边
switch (_case)
{
case 4:
case 10:
// 计算立方体的对角线上的点之间的线性插值
a = (_cube[4] - _cube[0]) * (_cube[6] - _cube[2]) - (_cube[7] - _cube[3]) * (_cube[5] - _cube[1]);
b = _cube[2] * (_cube[4] - _cube[0]) + _cube[0] * (_cube[6] - _cube[2])
- _cube[1] * (_cube[7] - _cube[3]) - _cube[3] * (_cube[5] - _cube[1]);
t = -b / (2 * a);
// 判断 t 是否在合法范围内
if (t < 0 || t > 1) return s > 0;
// 计算插值点的坐标
At = _cube[0] + (_cube[4] - _cube[0]) * t;
Bt = _cube[3] + (_cube[7] - _cube[3]) * t;
Ct = _cube[2] + (_cube[6] - _cube[2]) * t;
Dt = _cube[1] + (_cube[5] - _cube[1]) * t;
break;
case 6:
case 7:
case 12:
case 13:
// 根据 _case 设置参考边
switch (_case)
{
case 6: edge = test6[_config][2]; break;
case 7: edge = test7[_config][4]; break;
case 12: edge = test12[_config][3]; break;
case 13: edge = tiling13_5_1[_config][_subconfig][0]; break;
}
// 根据参考边计算插值点的坐标
switch (edge)
{
case 0:
t = _cube[0] / (_cube[0] - _cube[1]);
At = 0;
Bt = _cube[3] + (_cube[2] - _cube[3]) * t;
Ct = _cube[7] + (_cube[6] - _cube[7]) * t;
Dt = _cube[4] + (_cube[5] - _cube[4]) * t;
break;
case 1:
t = _cube[1] / (_cube[1] - _cube[2]);
At = 0;
Bt = _cube[0] + (_cube[3] - _cube[0]) * t;
Ct = _cube[4] + (_cube[7] - _cube[4]) * t;
Dt = _cube[5] + (_cube[6] - _cube[5]) * t;
break;
case 2:
t = _cube[2] / (_cube[2] - _cube[3]);
At = 0;
Bt = _cube[1] + (_cube[0] - _cube[1]) * t;
Ct = _cube[5] + (_cube[4] - _cube[5]) * t;
Dt = _cube[6] + (_cube[7] - _cube[6]) * t;
break;
case 3:
t = _cube[3] / (_cube[3] - _cube[0]);
At = 0;
Bt = _cube[2] + (_cube[1] - _cube[2]) * t;
Ct = _cube[6] + (_cube[5] - _cube[6]) * t;
Dt = _cube[7] + (_cube[4] - _cube[7]) * t;
break;
case 4:
t = _cube[4] / (_cube[4] - _cube[5]);
At = 0;
Bt = _cube[7] + (_cube[6] - _cube[7]) * t;
Ct = _cube[3] + (_cube[2] - _cube[3]) * t;
Dt = _cube[0] + (_cube[1] - _cube[0]) * t;
break;
case 5:
t = _cube[5] / (_cube[5] - _cube[6]);
At = 0;
Bt = _cube[4] + (_cube[7] - _cube[4]) * t;
Ct = _cube[0] + (_cube[3] - _cube[0]) * t;
Dt = _cube[1] + (_cube[2] - _cube[1]) * t;
break;
case 6:
t = _cube[6] / (_cube[6] - _cube[7]);
At = 0;
Bt = _cube[5] + (_cube[4] - _cube[5]) * t;
Ct = _cube[1] + (_cube[0] - _cube[1]) * t;
Dt = _cube[2] + (_cube[3] - _cube[2]) * t;
break;
case 7:
t = _cube[7] / (_cube[7] - _cube[4]);
At = 0;
Bt = _cube[6] + (_cube[5] - _cube[6]) * t;
Ct = _cube[2] + (_cube[1] - _cube[2]) * t;
Dt = _cube[3] + (_cube[0] - _cube[3]) * t;
break;
case 8:
t = _cube[0] / (_cube[0] - _cube[4]);
At = 0;
Bt = _cube[3] + (_cube[7] - _cube[3]) * t;
Ct = _cube[2] + (_cube[6] - _cube[2]) * t;
Dt = _cube[1] + (_cube[5] - _cube[1]) * t;
break;
case 9:
t = _cube[1] / (_cube[1] - _cube[5]);
At = 0;
Bt = _cube[0] + (_cube[4] - _cube[0]) * t;
Ct = _cube[3] + (_cube[7] - _cube[3]) * t;
Dt = _cube[2] + (_cube[6] - _cube[2]) * t;
break;
case 10:
t = _cube[2] / (_cube[2] - _cube[6]);
At = 0;
Bt = _cube[1] + (_cube[5] - _cube[1]) * t;
Ct = _cube[0] + (_cube[4] - _cube[0]) * t;
Dt = _cube[3] + (_cube[7] - _cube[3]) * t;
break;
case 11:
t = _cube[3] / (_cube[3] - _cube[7]);
At = 0;
Bt = _cube[2] + (_cube[6] - _cube[2]) * t;
Ct = _cube[1] + (_cube[5] - _cube[1]) * t;
Dt = _cube[0] + (_cube[4] - _cube[0]) * t;
break;
default:
printf("无效的边 %d\n", edge);
print_cube();
break;
}
break;
default:
printf("无效的模糊情况 %d\n", _case);
print_cube();
break;
}
// 根据计算结果判断立方体内部的状态
if (At >= 0) test++;
if (Bt >= 0) test += 2;
if (Ct >= 0) test += 4;
if (Dt >= 0) test += 8;
switch (test)
{
case 0: return s > 0;
case 1: return s > 0;
case 2: return s > 0;
case 3: return s > 0;
case 4: return s > 0;
case 5: if (At * Ct - Bt * Dt < FLT_EPSILON) return s > 0; break;
case 6: return s > 0;
case 7: return s < 0;
case 8: return s > 0;
case 9: return s > 0;
case 10: if (At * Ct - Bt * Dt >= FLT_EPSILON) return s > 0; break;
case 11: return s < 0;
case 12: return s > 0;
case 13: return s < 0;
case 14: return s < 0;
case 15: return s < 0;
}
return s < 0;
}
//_____________________________________________________________________________
看一下是如何添加三角形的，
trig为三角形数组，每三个边索引表示一个三角形如上面的图，nt为三角形个数
注意tv中是在init_temps中预计算的顶点索引（顶点数组在compute_intersection_points计算交点时就建好了）
void MarchingCubes::add_triangle(const char* trig, char n, int v12)
{
int tv[3]; // 存储三角形的三个顶点
for (int t = 0; t < 3 * n; t++)
{
switch (trig[t])
{
case 0: tv[t % 3] = get_x_vert(_i, _j, _k); break;
case 1: tv[t % 3] = get_y_vert(_i + 1, _j, _k); break;
case 2: tv[t % 3] = get_x_vert(_i, _j + 1, _k); break;
case 3: tv[t % 3] = get_y_vert(_i, _j, _k); break;
case 4: tv[t % 3] = get_x_vert(_i, _j, _k + 1); break;
case 5: tv[t % 3] = get_y_vert(_i + 1, _j, _k + 1); break;
case 6: tv[t % 3] = get_x_vert(_i, _j + 1, _k + 1); break;
case 7: tv[t % 3] = get_y_vert(_i, _j, _k + 1); break;
case 8: tv[t % 3] = get_z_vert(_i, _j, _k); break;
case 9: tv[t % 3] = get_z_vert(_i + 1, _j, _k); break;
case 10: tv[t % 3] = get_z_vert(_i + 1, _j + 1, _k); break;
case 11: tv[t % 3] = get_z_vert(_i, _j + 1, _k); break;
case 12: tv[t % 3] = v12; break;
default: break;
}
// 检查顶点是否合法
if (tv[t % 3] == -1)
{
printf("Marching Cubes: invalid triangle %d\n", _ntrigs + 1);
print_cube();
}
// 每三次处理一个三角形
if (t % 3 == 2)
{
// 如果当前三角形数组已满，则扩展数组 ***这个过程涉及到内存分配和拷贝，比较耗时
if (_ntrigs >= _Ntrigs)
{
Triangle *temp = _triangles;
_triangles = new Triangle[2 * _Ntrigs];
memcpy(_triangles, temp, _Ntrigs * sizeof(Triangle));
delete[] temp;
printf("%d allocated triangles\n", _Ntrigs);
_Ntrigs *= 2;
}
// 添加三角形到数组
Triangle *T = _triangles + _ntrigs++;
T->v1 = tv[0];
T->v2 = tv[1];
T->v3 = tv[2];
}
}
}
导出时，将顶点数组和三角形数组按相应格式写出就可以了
