本文已参与「新人创作礼」活动，一起开启掘金创作之路。

原文题目：Poisson surface reconstruction

点个star:github.com/Bailey-24/G…

一、实验目的

掌握用点云数据重建三角网格的计算方法

二、实验内容

实现泊松曲面重建算法泊松曲面重建属于隐函数方法实现。泊松表面重建的算法采取隐性拟合的方式，通过求解泊松方程来取得点云模型所描述的表面信息代表的隐性方程，通过对该方程进行等值面提取，从而得到具有几何实体信息的表面模型。

三、实验步骤

点个star:github.com/Bailey-24/G…

为了实现泊松曲面重建这个算法，我是按照总分结构去完成代码工作，也就是先补全 poisson surface reconstruction 文件，然后根据这个文件所需内容去完成剩下三个文件。

3.1 补全 poisson_surface_reconstruction.cpp 文件

3.1.1 功能需求

此函数输入点云和相应的法向向量，构建隐函数 g，然后用 marching cubes 算法输出三角网格的顶点和面。

3.1.2 步骤

a) 构建稀疏的三线性插值矩阵 Wx, Wy, Wz；

b) 通过 Wx, Wy, Wz 将给定的法线分布到 v 中的交错网格上，也就是 v 是被分配到交错网格上的法向向量上的拼接；

c) 构建和求解线性方程，得到隐函数g

d) 确定从 g 中提取的等值面；

e) 用 marching cubes 算法构建三角网格。

void poisson_surface_reconstruction( const Eigen::MatrixXd & P, const Eigen::MatrixXd & N, Eigen::MatrixXd & V, Eigen::MatrixXi & F) {
  ////////////////////////////////////////////////////////////////////////////
  // Construct FD grid, CONGRATULATIONS! You get this for free!
  ////////////////////////////////////////////////////////////////////////////

  // number of input points
  const int n = P.rows();
  // Grid dimensions
  int nx, ny, nz;
  // Maximum extent (side length of bounding box) of points
  double max_extent =
    (P.colwise().maxCoeff()-P.colwise().minCoeff()).maxCoeff();
  // padding: number of cells beyond bounding box of input points
  const double pad = 8;
  // choose grid spacing (h) so that shortest side gets 30+2*pad samples
  double h = max_extent/double(30+2*pad);
  // Place bottom-left-front corner of grid at minimum of points minus padding
  Eigen::RowVector3d corner = P.colwise().minCoeff().array()-pad*h;
  // Grid dimensions should be at least 3 
  nx = std::max((P.col(0).maxCoeff()-P.col(0).minCoeff()+(2.*pad)*h)/h,3.);
  ny = std::max((P.col(1).maxCoeff()-P.col(1).minCoeff()+(2.*pad)*h)/h,3.);
  nz = std::max((P.col(2).maxCoeff()-P.col(2).minCoeff()+(2.*pad)*h)/h,3.);
  // Compute positions of grid nodes
  Eigen::MatrixXd x(nx*ny*nz, 3);
  for(int i = 0; i < nx; i++) 
  {
    for(int j = 0; j < ny; j++)
    {
      for(int k = 0; k < nz; k++)
      {
         // Convert subscript to index
         const auto ind = i + nx*(j + k * ny);
         x.row(ind) = corner + h*Eigen::RowVector3d(i,j,k);
      }
    }
  }
  Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);

  ////////////////////////////////////////////////////////////////////////////
    Eigen::SparseMatrix<double> Wx, Wy, Wz;
    // staggered interpolation grids
    Eigen::RowVector3d corner_x = corner + Eigen::RowVector3d(h/2, 0,   0);
    Eigen::RowVector3d corner_y = corner + Eigen::RowVector3d(0,   h/2, 0);
    Eigen::RowVector3d corner_z = corner + Eigen::RowVector3d(0,   0,   h/2);
    fd_interpolate(nx-1, ny, nz, h, corner_x, P, Wx);
    fd_interpolate(nx, ny-1, nz, h, corner_y, P, Wy);
    fd_interpolate(nx, ny, nz-1, h, corner_z, P, Wz);
    
    // blend normals
    Eigen::VectorXd vx(nx*ny*nz - ny*nz), vy(nx*ny*nz - nx*nz), vz(nx*ny*nz - nx*ny);
    Eigen::VectorXd vxy(2*nx*ny*nz - ny*nz - nx*nz), v(3*nx*ny*nz - nx*ny - ny*nz - nx*nz);
    vx = Wx.transpose() * N.col(0);
    vy = Wy.transpose() * N.col(1);
    vz = Wz.transpose() * N.col(2);
    igl::cat(1, vx, vy, vxy);
    igl::cat(1, vxy, vz, v);

    // get gradient
    Eigen::SparseMatrix<double> G;
    fd_grad(nx, ny, nz, h, G);

    // G^T * G * g = G^T *v
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver(G.transpose() * G);
    g = solver.solve(G.transpose() * v);

    // determine sigma
    Eigen::SparseMatrix<double> W;
    fd_interpolate(nx, ny, nz, h, corner, P, W);
    double sigma = (W * g).sum()/n;
    g = g.array() - sigma;

  // Run black box algorithm to compute mesh from implicit function: this
  // function always extracts g=0, so "pre-shift" your g values by -sigma
  ////////////////////////////////////////////////////////////////////////////
  igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F);
}

3.2 补全 fd_interpolate 文件

3.2.1 功能需求

函数输入有限差分网格，它由网格的每一侧的节点数（nx, ny和nz），网格间距，处于最下最左最前的一个拐角坐标点和点云组成，输出一个三线性插值的点的权重矩阵 W 。

3.2.2 步骤

a) 将点的坐标归一化到拐角； b) 计算这个点的八个邻点的权重，其中 weight_list 是一个三元组（行，列，值），列的索引公式是i + j*n_x + k* n_y * n_x，值是三线性插值公式的系数，计算公式见参考文献1.

void fd_interpolate( const int nx, const int ny, const int nz, const double h, const Eigen::RowVector3d & corner, const Eigen::MatrixXd & P, Eigen::SparseMatrix<double> & W) {
	using namespace std;
	typedef Eigen::Triplet<double> T;
	vector<T> weight_list;
	weight_list.reserve(P.rows()*8*2);
	for (int row = 0; row < P.rows(); row++) {
		double dx = P(row, 0) - corner(0);
		double dy = P(row, 1) - corner(1);
		double dz = P(row, 2) - corner(2);
		int i = int(floor(dx / h));
		double x_value = fmod(dx, h) / h;
		int j = int(floor(dy / h));
		double y_value = fmod(dy, h) / h;
		int k = int(floor(dz / h));
		double z_value = fmod(dz, h) / h;

		// compute the value for the 8 neighboring nodes
		weight_list.push_back(T(row, i + nx * (j + k * ny), 
			(1.0 - x_value) * (1.0 - y_value) * (1.0 - z_value)));
		weight_list.push_back(T(row, (i + 1) + nx * (j + k * ny),
			x_value * (1.0 - y_value) * (1.0 - z_value)));
		weight_list.push_back(T(row, i + nx * ((j + 1) + k * ny),
			(1.0 - x_value) * y_value * (1.0 - z_value)));
		weight_list.push_back(T(row, (i + 1) + nx * ((j + 1) + k * ny),
			x_value * y_value * (1.0 - z_value)));
		weight_list.push_back(T(row, i + nx * (j + (k + 1) * ny),
			(1.0 - x_value) * (1.0 - y_value) * z_value));
		weight_list.push_back(T(row, (i + 1) + nx * (j + (k + 1) * ny),
			x_value * (1.0 - y_value) * z_value));
		weight_list.push_back(T(row, i + nx * ((j + 1) + (k + 1) * ny),
			(1.0 - x_value) * y_value * z_value));
		weight_list.push_back(T(row, (i + 1) + nx * ((j + 1) + (k + 1) * ny),
			x_value * y_value * z_value));
	}
	W.resize(P.rows(), nx*ny*nz);
	W.setFromTriplets(weight_list.begin(), weight_list.end());
}

3.3 补全 fd_grad 文件

3.3.1 功能需求

函数输入有限差分网格，它由网格的每一侧的节点数（nx, ny和nz）和网格间距组成，计算每个分量在其各自的交错网格上的梯度，将其构造成稀疏矩阵 G 并输出。

3.3.2 步骤

首先通过下面的函数计算出每个方向上的梯度，然后将其拼接在一起，构成稀疏矩阵 G 。

void fd_grad( const int nx, const int ny, const int nz, const double h, Eigen::SparseMatrix<double> & G) {
  Eigen::SparseMatrix<double> Dx, Dy, Dz;
  fd_partial_derivative(nx, ny, nz,h, 0, Dx);
  fd_partial_derivative(nx, ny, nz,h, 1, Dy);
  fd_partial_derivative(nx, ny, nz,h, 2, Dz);

  // Attempted to use coeffRef to populate G from ground up but I was again stuck on assertion failure
  // Reference Sarah's approach on using igl::cat
  Eigen::SparseMatrix<double> Dxy((nx - 1) * ny * nz + nx * (ny - 1) * nz, nx * ny * nz);
  G.resize(Dxy.size() + nx * ny * (nz - 1), nx * ny * nz);
  igl::cat(1, Dx, Dy, Dxy);
  igl::cat(1, Dxy, Dz, G);
}

3.4 补全 fd_partial_derivative 文件

3.4.1 功能需求

函数输入有限差分网格，它由网格的每一侧的节点数（nx, ny和nz），网格间距和指定的方向，输出稀疏矩阵 D，它是在指定方向上的交错网格的一阶偏导数。

3.4.2 步骤

a) 找到在交错网格上每个点； b) 对这个点的左边的第一个常规网格赋值 -1； c) 对这个点的右边的第一个常规网格赋值 1； d) 构造稀疏矩阵 D。

void fd_partial_derivative( const int nx, const int ny, const int nz, const double h, const int dir, Eigen::SparseMatrix<double> & D) {
    int nx_D = nx, ny_D = ny, nz_D = nz;

    switch (dir)
    {
    case 0:
        nx_D--;
        break;
    case 1:
        ny_D--;
        break;
    case 2:
        nz_D--;
        break;

    default:
        assert(0);
    }
    // resize D so that each row corresponds to a node in the staggered grid
    // and each column corresponds to a node in the regular grid
    D.resize(nx_D * ny_D * nz_D, nx * ny * nz);

    // loop through each of the nodes in the staggered grid 
    // and assign -1 to the regular grid node "left"-neighbour
    // and 1 to the regular grid node "right"-neighbour
    for (int i=0; i<nx_D; i++)
    {
        for (int j=0; j<ny_D; j++)
        {
            for (int k=0; k<nz_D; k++)
            {
                int staggered_ijk = i + j * nx_D + k * ny_D * nx_D;
                int l_prev = i + j * nx + k * ny * nx;
                int l_curr;

                switch (dir)
                {
                case 0:
                    l_curr = (i + 1) + j * nx + k * ny * nx;
                    break;
                case 1:
                    l_curr = i + (j + 1) * nx + k * ny * nx;
                    break;
                case 2:
                    l_curr = i + j * nx + (k + 1) * ny * nx;
                    break;
                }
                D.insert(staggered_ijk, l_prev) = -1;
                D.insert(staggered_ijk, l_curr) = 1;
            }
        }
    }
}

3.4 结果

计算出的结果与作者相同，视觉效果也不错。

四、总结分析

坐标系很关键在补全三线性插值函数时，如果坐标系不对，则计算的系数也会不对。
索引很重要在补全三线性插值函数和偏导数函数时，常常因为索引不对，而得不到正确的结果。

五、参考文献

【1】线性插值与双/三线性插值

今天的文章【传统图形学算法】泊松曲面重建将点云重建成表面 | 传统重建算法的佼佼者分享到此就结束了，感谢您的阅读。

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容，请发送邮件至举报，一经查实，本站将立刻删除。
如需转载请保留出处：https://bianchenghao.cn/19912.html

【传统图形学算法】泊松曲面重建将点云重建成表面 | 传统重建算法的佼佼者