PCL源码剖析之MarchingCubes算法

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;
  }
}

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 ());
      }
}

点q处的法向量是单位法向量，所以点q到切平面的距离是dot(N(p), vec(p, q))

从代码中可以看出，这里的K = 1，即求出最近邻点。

文章参考于：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

PCL源码剖析之MarchingCubes算法

原文：http://blog.csdn.net/lming_08/article/details/19432877