本文已参加「新人创造礼」活动， 一起开启掘金创造之路。
原文题目：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) 构建和求解线性方程
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】 线性插值与双/三线性插值
