MarchingCubes算法简介
MarchingCubes(移动立方体)算法是目前三围数据场等值面生成中最常用的方法。它实际上是一个分而治之的方法,把等值面的抽取分布于每个体素中进行。对于每个被处理的体素,以三角面片逼近其内部的等值面片。每个体素是一个小立方体,构造三角面片的处理过程对每个体素都“扫描”一遍,就好像一个处理器在这些体素上移动一样,由此得名移动立方体算法。
MC算法主要有三步:1.将点云数据转换为体素网格数据;2.使用线性插值对每个体素抽取等值面;3.对等值面进行网格三角化
PCL源码剖析之MarchingCubesHoppe
PCL中使用MarchingCubesHoppe类进行三维重建执行的函数体为performReconstruction(),其代码如下:
-
template <
typename PointNT>
void
-
pcl::MarchingCubes<PointNT>::performReconstruction (pcl::PolygonMesh &output)
-
{
-
if (!(iso_level_ >=
0 && iso_level_ <
1))
-
{
-
PCL_ERROR (
"[pcl::%s::performReconstruction] Invalid iso level %f! Please use a number between 0 and 1.\n", getClassName ().c_str (), iso_level_);
-
output.cloud.width = output.cloud.height =
0;
-
output.cloud.data.clear ();
-
output.polygons.clear ();
-
return;
-
}
-
-
// Create grid
-
grid_ =
std::
vector<
float> (res_x_*res_y_*res_z_,
0.0f);
-
-
// Populate tree
-
tree_->setInputCloud (input_);
-
-
getBoundingBox ();
-
-
// Transform the point cloud into a voxel grid
-
// This needs to be implemented in a child class
-
voxelizeData ();
-
-
-
-
// Run the actual marching cubes algorithm, store it into a point cloud,
-
// and copy the point cloud + connectivity into output
-
pcl::PointCloud<PointNT> cloud;
-
-
for (
int x =
1; x < res_x_
-1; ++x)
-
for (
int y =
1; y < res_y_
-1; ++y)
-
for (
int z =
1; z < res_z_
-1; ++z)
-
{
-
Eigen::
Vector3i index_3d (x, y, z);
-
std::
vector<
float> leaf_node;
-
getNeighborList1D (leaf_node, index_3d);
-
createSurface (leaf_node, index_3d, cloud);
-
}
-
pcl::toPCLPointCloud2 (cloud, output.cloud);
-
-
output.polygons.resize (cloud.size () /
3);
-
for (
size_t i =
0; i < output.polygons.size (); ++i)
-
{
-
pcl::Vertices v;
-
v.vertices.resize (
3);
-
for (
int j =
0; j <
3; ++j)
-
v.vertices[j] =
static_cast<
int> (i) *
3 + j;
-
output.polygons[i] = v;
-
}
-
}
可以看出PCL将会产生res_x_ * res_y_ * res_z_个网格,即为Resolution分辨率。voxelizeData ();即将点云数据转换为体素网格数据,其实现如下:
-
template <
typename PointNT>
void
-
pcl::MarchingCubesHoppe<PointNT>::voxelizeData ()
-
{
-
for (
int x =
0; x < res_x_; ++x)
-
for (
int y =
0; y < res_y_; ++y)
-
for (
int z =
0; z < res_z_; ++z)
-
{
-
std::
vector<
int> nn_indices;
-
std::
vector<
float> nn_sqr_dists;
-
-
Eigen::Vector3f point;
-
point[
0] = min_p_[
0] + (max_p_[
0] - min_p_[
0]) *
float (x) /
float (res_x_);
-
point[
1] = min_p_[
1] + (max_p_[
1] - min_p_[
1]) *
float (y) /
float (res_y_);
-
point[
2] = min_p_[
2] + (max_p_[
2] - min_p_[
2]) *
float (z) /
float (res_z_);
-
-
PointNT p;
-
p.getVector3fMap () = point;
-
-
tree_->nearestKSearch (p,
1, nn_indices, nn_sqr_dists);
-
-
grid_[x * res_y_*res_z_ + y * res_z_ + z] = input_->points[nn_indices[
0]].getNormalVector3fMap ().dot (
-
point - input_->points[nn_indices[
0]].getVector3fMap ());
-
}
-
}
该函数对每个体素网格数据进行赋值,其值为一符号距离函数值,其定义为:f(Pi) = (Pi - Oi) * N(Oi), 这里Pi为给定的点,Oi为Pi周围K近邻点集(输入点云的子集)的中心, N(Oi)为点Oi的法向量,中间的*为数量积;求出的值其实是点Oi到过点Pi的有向切平面的距离,图示如下:
点q处的法向量是单位法向量,所以点q到切平面的距离是dot(N(p), vec(p, q))
从代码中可以看出,这里的K = 1,即求出最近邻点。
下面的代码描述了对每个体素网格的处理过程,其主要过程是计算出每个体素网格与等值面的交点,然后按一定顺序将交点连接,从而形成三角面片。
-
for (
int x =
1; x < res_x_
-1; ++x)
-
for (
int y =
1; y < res_y_
-1; ++y)
-
for (
int z =
1; z < res_z_
-1; ++z)
-
{
-
Eigen::
Vector3i index_3d (x, y, z);
-
std::
vector<
float> leaf_node;
-
getNeighborList1D (leaf_node, index_3d);
-
createSurface (leaf_node, index_3d, cloud);
-
}
getNeighorList1D(leaf_node, index_3d);即是求出当前体素网格的8个顶点对应符号距离函数值,即数组grid_中对应的值。实现代码如下:
-
template <
typename PointNT>
void
-
pcl::MarchingCubes<PointNT>::getNeighborList1D (
std::
vector<
float> &leaf,
-
Eigen::Vector3i &index3d)
-
{
-
leaf =
std::
vector<
float> (
8,
0.0f);
-
-
leaf[
0] = getGridValue (index3d);
-
leaf[
1] = getGridValue (index3d + Eigen::Vector3i (
1,
0,
0));
-
leaf[
2] = getGridValue (index3d + Eigen::Vector3i (
1,
0,
1));
-
leaf[
3] = getGridValue (index3d + Eigen::Vector3i (
0,
0,
1));
-
leaf[
4] = getGridValue (index3d + Eigen::Vector3i (
0,
1,
0));
-
leaf[
5] = getGridValue (index3d + Eigen::Vector3i (
1,
1,
0));
-
leaf[
6] = getGridValue (index3d + Eigen::Vector3i (
1,
1,
1));
-
leaf[
7] = getGridValue (index3d + Eigen::Vector3i (
0,
1,
1));
-
}
createSurface (leaf_node, index_3d, cloud);即是求出每个体素网格与等值面的交点,然后按一定顺序将交点连接,从而形成三角面片。下面分片段剖析createSurface ()函数代码。
-
int cubeindex =
0;
-
Eigen::Vector3f vertex_list[
12];
-
if (leaf_node[
0] < iso_level_) cubeindex |=
1;
-
if (leaf_node[
1] < iso_level_) cubeindex |=
2;
-
if (leaf_node[
2] < iso_level_) cubeindex |=
4;
-
if (leaf_node[
3] < iso_level_) cubeindex |=
8;
-
if (leaf_node[
4] < iso_level_) cubeindex |=
16;
-
if (leaf_node[
5] < iso_level_) cubeindex |=
32;
-
if (leaf_node[
6] < iso_level_) cubeindex |=
64;
-
if (leaf_node[
7] < iso_level_) cubeindex |=
128;
此段代码是将8个顶点的标量值与等值面相比较,如果标量值小于等值面值(即顶点在等值面下面),则将cubeindex相应的位置为1。这样就可以知道8个顶点中哪些在等值面下,哪些在等值面之上了。
立方体中顶点与棱边的编号如下所示:
例如,如果顶点3的值在等值面值之下并且所有其他顶点的值都在等值面之上,那么我们可以通过切割边2、3、11来创建一个三角面片。
算法使用一个边表将cubeindex映射为一个12bit的数值,每一位与一条边相关,如果边没有被等值面切割则设为0,切割则设为1。如果没有边被切割那么表返回0,这种情况发生在当cubeindex = 0(所有顶点在等值面之下)或0xff(所有顶点在等值面之上)。举个之前的例子,如果只有顶点3在等值面下面,cubeindex将会等于0000 1000 或8。边表edgeTable定义在marching_cubes.h文件中:
-
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
-
};
edgeTable[8] = 0x80c = 1000 0000 1100。这就表示边2、3、11与等值面相交。
判断出边与等值面相交之后,就要确定具体的交点。这里PCL使用线性插值,线性插值具体见线性插值。代码如下:
-
// Find the vertices where the surface intersects the cube
-
if (edgeTable[cubeindex] &
1)
-
interpolateEdge (p[
0], p[
1], leaf_node[
0], leaf_node[
1], vertex_list[
0]);
-
if (edgeTable[cubeindex] &
2)
-
interpolateEdge (p[
1], p[
2], leaf_node[
1], leaf_node[
2], vertex_list[
1]);
-
if (edgeTable[cubeindex] &
4)
-
interpolateEdge (p[
2], p[
3], leaf_node[
2], leaf_node[
3], vertex_list[
2]);
-
if (edgeTable[cubeindex] &
8)
-
interpolateEdge (p[
3], p[
0], leaf_node[
3], leaf_node[
0], vertex_list[
3]);
-
if (edgeTable[cubeindex] &
16)
-
interpolateEdge (p[
4], p[
5], leaf_node[
4], leaf_node[
5], vertex_list[
4]);
-
if (edgeTable[cubeindex] &
32)
-
interpolateEdge (p[
5], p[
6], leaf_node[
5], leaf_node[
6], vertex_list[
5]);
-
if (edgeTable[cubeindex] &
64)
-
interpolateEdge (p[
6], p[
7], leaf_node[
6], leaf_node[
7], vertex_list[
6]);
-
if (edgeTable[cubeindex] &
128)
-
interpolateEdge (p[
7], p[
4], leaf_node[
7], leaf_node[
4], vertex_list[
7]);
-
if (edgeTable[cubeindex] &
256)
-
interpolateEdge (p[
0], p[
4], leaf_node[
0], leaf_node[
4], vertex_list[
8]);
-
if (edgeTable[cubeindex] &
512)
-
interpolateEdge (p[
1], p[
5], leaf_node[
1], leaf_node[
5], vertex_list[
9]);
-
if (edgeTable[cubeindex] &
1024)
-
interpolateEdge (p[
2], p[
6], leaf_node[
2], leaf_node[
6], vertex_list[
10]);
-
if (edgeTable[cubeindex] &
2048)
-
interpolateEdge (p[
3], p[
7], leaf_node[
3], leaf_node[
7], vertex_list[
11]);
文章参考于:http://www.doc88.com/p-5475997688638.html
http://paulbourke.net/geometry/polygonise/
http://books.google.com.hk/books?id=4k4kvDwP-lgC&printsec=frontcover&hl=zh-CN#v=onepage&q&f=false