MC（移动立方体）算法

LV小猪精

已于 2024-02-22 15:44:01 修改

阅读量5.3k

点赞数 5

分类专栏：算法计算机图形学文章标签： MC算法移动立方体网格重构

于 2022-03-23 22:11:26 首次发布

本文链接：https://blog.csdn.net/weixin_44478077/article/details/123696173

版权

Marching Cubes（移动立方体）

Marching squares 算法
Marching Cubes算法

📘写在最前面：
在学习MC算法之前先看看二维的MS算法，MC算法是二维MS算法的拓展，通过学习MS的基本思路能够更快学习MC算法，因此在这里我们首先提一下MS的基本算法流程。

Marching squares 算法

MS算法的主要流程：

  1. 遍历每个栅格，从16种固定情况中选择类别

  2. 利用线性插值，结合网格点数值找寻交点

  3. 根据线性插值的结果连线

具体的操作过程如图所示：
在这里插入图片描述

Marching Cubes算法

MC算法也被称作“等值面提取（Isosurface Extraction）”是面绘制算法中的经典算法，算法的主要精髓为：在三维离散数据场中通过线性差值来逼近等值面。
在学习MC算法前需要给定一些基本概念：首先定义一个立方体单元为一个体素，每一个体元均由8个顶点所构成。

体素顶点有两种不同的状态量所表示：

高于或等于势值表示在物体表面的内部
低于势值表示在物体表面的外部
因此，体素的一个顶点有两中可能的状态，则一个体素（8个顶点）就一共有2^8=256种状态。根据旋转、映射不变等特性，体素的状态可以被归纳为以下15种基本构型：

在这里插入图片描述
为了存储方便，将其记录在表中（如下所示）

int MarchingCubes::ms_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   };

对模型包围盒空间划分体,从下到上每个体素单元均遍历一次
在这里插入图片描述
遍历完成后得到构建的模型

构建得到的三角网格模型如图所示

附上c++实现的完整代码：

// utils.h
#pragma once
#include <math.h>

#define NULL 0
#define BYTE unsigned char

template <class T> class VECTOR3 {
   
public:
	T x;
	T y;
	T z;

	VECTOR3(): x(0), y(0),z(0) {
   }
	VECTOR3(const VECTOR3<int> &other): x(other.x), y(other.y), z(other.z) {
   }
	VECTOR3(const T _x, const T _y, const T _z) : x(_x), y(_y), z(_z) {
   }

	bool operator == ( const VECTOR3<T>& other ) const {
   return (other.x==x && other.y==y && other.z==z); }
	bool operator != ( const VECTOR3<T>& other ) const {
   return (other.x!=x || other.y!=y || other.z!=z); } 

    // binary operators with scalars
	VECTOR3<T> operator + ( T scalar ) const {
   return VECTOR3<T>(x+scalar,y+scalar,z+scalar);}
	VECTOR3<T> operator - ( T scalar ) const {
   return VECTOR3<T>(x-scalar,y-scalar,z-scalar);}
	VECTOR3<T> operator * ( T scalar ) const {
   return VECTOR3<T>(x*scalar,y*scalar,z*scalar);}
	VECTOR3<T> operator / ( T scalar ) const {
   return VECTOR3<T>(x/scalar,y/scalar,z/scalar);}

	// binaray operpators  with vectors
	VECTOR3<T> operator + ( const VECTOR3<T>& other ) const {
   return VECTOR3<T>(x+other.x,y+other.y,z+other.z);}
	VECTOR3<T> operator - ( const VECTOR3<T>& other ) const {
   return VECTOR3<T>(x-other.x,y-other.y,z-other.z);}
	VECTOR3<T> operator * ( const VECTOR3<T>& other ) const {
   return VECTOR3<T>(x*other.x,y*other.y,z*other.z);}
	VECTOR3<T> operator / ( const VECTOR3<T>& other ) const {
   return VECTOR3<T>(x/other.x,y/other.y,z/other.z);}

	T maxVal() {
   return (x>y) ? ((x>z) ? x : z) : ((y>z) ? y : z);}
	T minVal() {
   return (x<y) ? ((x<z) ? x : z) : ((y<z) ? y : z);}
	T volume() {
   return x*y*z;}
	T length() {
   return sqrt(T(x*x+y*y+z*z));}
	void normalize() {
   T len = length(); x/=len;y/=len;;z/=len;}
	void normalize(float epsilon) {
   
		T len = length();
		if (len > epsilon) {
   
			x/=len;
			y/=len;
			z/=len;
		} else {
    // specify some arbitrary normal
			x = 0;
			y = 0;
			z = 1;
		}
	}
};


void loadData(char *fileName, bool bIsFloat, VECTOR3<int>* iVolSize, float **data);
void loadData(char *fileName, bool bIsFloat, VECTOR3<int> iVolSize, float **data);

//utils.cpp
#include "utils.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>


void loadData(char *fileName, bool bIsFloat, VECTOR3<int>* iVolSize, float **data) {
   
	printf("\nTrying to load %s Dataset...\n",fileName);

	FILE *stream;	
	if((stream = fopen(fileName, "rb"))==NULL) {
   
		printf("Error while loading data from file %s.\n",fileName);
		exit(0);
	} else {
   
		fread(iVolSize , sizeof(VECTOR3<int>), 1, stream);
		*data = new float[iVolSize->volume()];
		if (bIsFloat) {
   
			fread(*data, sizeof(float), iVolSize->volume(), stream);
		} else {
   
			BYTE* temp_data = new BYTE[iVolSize->volume()];
			fread(temp_data, sizeof(BYTE), iVolSize->volume(), stream);
			for (int i = 0;i<iVolSize->volume();i++) (*data)[i] = float(temp_data[i])/255.0f;
			delete [] temp_data;
		}

		fclose( stream );
		printf("\nData loaded\n");
	}
}

void loadData(char *fileName, bool bIsFloat, VECTOR3<int> iVolSize, float **data) {
   
	printf("\nTrying to load %s Dataset...\n",fileName);

	FILE *stream;	
	if((stream = fopen(fileName, "rb"))==NULL) {
   
		printf("Error while loading data from file %s.\n",fileName);
		exit(0);
	} else {
   
		*data = new float[iVolSize.volume()];
		if (bIsFloat) {
   
			fread(*data, sizeof(float), iVolSize.volume(), stream);
		} else {
   
			BYTE* temp_data = new BYTE[iVolSize.volume()];
			fread(temp_data, sizeof(BYTE), iVolSize.volume(), stream);
			for (int i = 0;i<iVolSize.volume();i++) (*data)[i] = float(temp_data[i])/255.0f;
			delete [] temp_data;
		}

		fclose( stream );
		printf("\nData loaded\n");
	}
}


//MarchingCube.h
#pragma once

#include <memory.h>
#include "utils.h"

#define EPSILON 0.000001f
#define DATA_INDEX(I, J, K, IDIM, JDIM) ((I) + ((IDIM) * (J)) + ((IDIM * JDIM) * (K)))
#define EDGE_INDEX(E, I, J, IDIM) ((E) + (12 * (I)) + ((12 * (IDIM)) * (J)))
#define NO_EDGE -1

class LayerTempData 
{
   
private:
	VECTOR3<int> m_vVolSize;
public:
	LayerTempData(VECTOR3<int> vVolSize, float* pfVolume) 
	{
   
		m_vVolSize = vVolSize;
		
		pfBotData	= pfVolume;
		pfTopData	= pfVolume+DATA_INDEX(0, 0, 1, vVolSize.x, vVolSize.y);
		piEdges		= new int[(vVolSize.x-1) * (vVolSize.y-1) * 12];  // allocate storage to hold the indexing tags for edges in the layer
		
		for (int i = 0; i < (vVolSize.x-1) * (vVolSize.y-1) * 12; i++)	piEdges[i] = NO_EDGE;  // init edge list
	}
	
	virtual ~LayerTempData() 
	{
   
		delete [] piEdges;
		delete [] pfTopData;
		delete [] pfBotData;
	}
	
	void nextIteration() 
	{
   
		// update the layer for this iteration
		pfBotData = pfTopData;
		// now topData points to next layer of scalar data
		pfTopData += DATA_INDEX(0, 0, 1, m_vVolSize.x, m_vVolSize.y);
		// percolate the last layer's top edges to this layer's bottom edges
		for (int iY = 0; iY < m_vVolSize.y-1; iY++) 
		{
   
			for (int iX = 0; iX < m_vVolSize.x-1; iX++) 
			{
   
				for (int iEdges = 0; iEdges < 4; iEdges++) 
				{
   
					piEdges[EDGE_INDEX(iEdges, iX, iY, m_vVolSize.x-1)] =  piEdges[EDGE_INDEX(iEdges+4, iX, iY, m_vVolSize.x-1)];
				}
				// reinitialize all of the remaining edges
				for (int  iEdges = 4; iEdges < 12; iEdges++) 
				{
   
					piEdges[EDGE_INDEX(iEdges, iX, iY, m_vVolSize.x-1)] = NO_EDGE;
				}
			}
		}
	}
	
	float*	pfBotData;
	float*	pfTopData;
	int*	piEdges;  // tag indexing into vertex list
};

class Isosurface{
   
public:
	VECTOR3<float>*	vfVertices;
	VECTOR3<float>*	vfNormals;
	VECTOR3<int>*	viTriangles;
	int				iVertices;
	int				iTriangles;

	Isosurface() {
   
	    vfVertices  = NULL;
	    vfNormals   = NULL;
	    viTriangles = NULL;
	    iVertices   = 0;
	    iTriangles  = 0;
	}

	Isosurface(int iMaxVertices, int iMaxTris) 
	{
   
	    vfVertices  = new VECTOR3<float>[iMaxVertices];
	    vfNormals   = new VECTOR3<float>[iMaxVertices];
	    viTriangles = new VECTOR3<int>[iMaxTris];
	    iVertices   = 0;
	    iTriangles  = 0;
	}

	virtual ~Isosurface() 
	{
   
		delete [] vfVertices;
		delete [] vfNormals;
		delete [] viTriangles;
	}

	int AddTriangle(int a, int b, int c) 
	{
   
		viTriangles[iTriangles++] = VECTOR3<int>(a,b,c);
		return iTriangles-1;
	}

	int AddVertex(VECTOR3<float> v, VECTOR3<float> n) 
	{
   
		vfVertices[iVertices] = v;
		vfNormals[iVertices++] = n; 
		return iVertices-1;
	}

	void AppendData(const Isosurface* other) 
	{
   
		// if verts in other, expand the storage this surface
		if (other->iVertices > 0) 
		{
   

			// create new mem
			VECTOR3<float>*	temp_Vertices  = new VECTOR3<float>[iVertices  + other->iVertices];
			VECTOR3<float>*	temp_Normals   = new VECTOR3<float>[iVertices  + other->iVertices];
			VECTOR3<int>*	temp_Triangles = new VECTOR3<int>[iTriangles + other->iTriangles];

			// copy "old" data
			memcpy(temp_Vertices, vfVertices, sizeof(VECTOR3<float>)*iVertices);
			memcpy(temp_Normals, vfNormals, sizeof(VECTOR3<float>)*iVertices);
			memcpy(temp_Triangles, viTriangles, sizeof(VECTOR3<int>)*iTriangles);

			// append new data
			memcpy(temp_Vertices+iVertices, other->vfVertices, sizeof(VECTOR3<float>)*other->iVertices);
			memcpy(temp_Normals+iVertices, other->vfNormals, sizeof(VECTOR3<float>)*other->iVertices);
			memcpy(temp_Triangles+iTriangles, other->viTriangles, sizeof(VECTOR3<int>)*other->iTriangles);

			// delete "old" data
			delete [] vfVertices;
			delete [] vfNormals;
			delete [] viTriangles;

			// rename
			vfVertices  = temp_Vertices;
			vfNormals   = temp_Normals;
			viTriangles = temp_Triangles;
		}

		// update this list's counters
		iVertices  += other->iVertices;
		iTriangles += other->iTriangles;
	}
};


class MarchingCubes
{
   
protected:
	static int		ms_edgeTable[256];
	static int		ms_triTable[256][16];

	float			m_fFlipFact;
	float			m_bFlipSurface;
	VECTOR3<int>	m_vVolSize;
	float*			m_pfVolume;
	float			m_fIsoValue;

	virtual void MarchLayer(LayerTempData *layer, int iLayer);
	virtual int MakeVertex(int whichEdge, int i, int j, int k, Isosurface* sliceIso);
	virtual VECTOR3<float> InterpolateNormal(float fValueAtPos, VECTOR3<int> vPosition);

public:
	Isosurface*		m_Isosurface;
	VECTOR3<float>  BoxMin;     //网格顶点中最左下点坐标
	float           spacing;    //网格间距，这里假定三个方向的间距一致

	MarchingCubes(void);
	virtual ~MarchingCubes(void);

	virtual void SetVolume(int iSizeX, int iSizeY, int iSizeZ, float* pfVolume);
	virtual void Process(float fIsoValue, bool bFlipSurface=false);
};

//MarchingCube.cpp
#include "marchingcubes.h"
int MarchingCubes::ms_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

最低0.47元/天解锁文章

LV小猪精

关注

5
点赞
踩
53

收藏

觉得还不错? 一键收藏
打赏
7
评论
MC（移动立方体）算法

Marching Cubes（移动立方体）MC算法MC算法写在最前面 MS算法MC算法也被称作“等值面提取（Isosurface Extraction）”是面绘制算法中的经典算法，算法的主要精髓为：在三维离散数据场中通过线性差值来逼近等值面。在学习MC算法前需要给定一些基本概念：首先定义一个立方体单元为一个体素，每一个体元均由8个顶点所构成。体素顶点有两种不同的状态量所表示：1.高于或等于势值表示在物体表面的内部2.低于势值表示在物体表面的外部因此，体素的一个顶点有两中可能的状态，则一个
复制链接

扫一扫