Marching Cubes 算法再探

本文主要是介绍Marching Cubes 算法再探，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

之前做NeRF相关工作时简单看过，但没有深究其实现，知其然不知其所以然的程度，算是初探。

基础

为了对MC有大致了解，可以先根据Marching Cubes 算法初探，理解一下Marching Cubes这个名字的由来，其二维空间中的示例可谓一目了然。

接下来结合Efficient Implementation of Marching Cubes’ Cases with Topological Guarantees
的论文和代码，来探究一下CPU版本经典MC的具体实现过程。

过程

整个过程的核心部分可以抽象为以下步骤：

1. 对每一个cube，比较其顶点与等值面的大小，得到顶点的8bit案例编码

顶点1、5、7值为正，则得到10100010，即162

2. 根据顶点案例的8bit编码在查找表中找到相应三角形数组

162对应的三角形数组为{5,0,1,5,4,0,7,6,11}

4. 根据三角形数组添加\绘制三角形

这个三角形数组每三个值代表一个三角形的三个顶点位置（先以cube边索引形式给出）

添加三角形[5,0,1]

添加三角形[5,4,0]

添加三角形[7,6,11]

4. 遍历所有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 配置elseadd_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 配置elseadd_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.2else{if( test_interior( test6[_config][1]) )add_triangle( tiling6_1_1[_config], 3 ) ; // 6.1.1else{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 ) ;elseadd_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.1else{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.1elseadd_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.1else{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.1elseadd_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];}}
}

导出时，将顶点数组和三角形数组按相应格式写出就可以了

---

这篇关于Marching Cubes 算法再探的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---