diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..be48f30 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,13 +1,39 @@ #include "fd_grad.h" -void fd_grad( - const int nx, - const int ny, - const int nz, - const double h, - Eigen::SparseMatrix & G) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// +typedef Eigen::Triplet T; + +void intervalue2(std::vector &coef, int row, int aidx, int bidx) { + T tmp(row, aidx, -1); + T tmp2(row, bidx, 1); + coef.push_back(tmp); + coef.push_back(tmp2); } + +extern int gidx(double i, double j, double k, int nx, int ny, int nz); + +void fd_grad(const int nx, const int ny, const int nz, const double h, + Eigen::SparseMatrix & G) { + //////////////////////////////////////////////////////////////////////////// + // Add your code here + + // first order + std::vector coef; + for (int i = 0; i < nx - 1; i++) + for (int j = 0; j < ny - 1; j++) + for (int k = 0; k < nz - 1; k++) { + + int st = gidx(i, j, k, nx, ny, nz); + int edx = gidx(i + 1, j, k, nx, ny, nz); + int edy = gidx(i, j + 1, k, nx, ny, nz); + int edz = gidx(i, j, k + 1, nx, ny, nz); + + int row = st * 3; + intervalue2(coef, row, st, edx); + intervalue2(coef, row + 1, st, edy); + intervalue2(coef, row + 2, st, edz); + } + + G.setFromTriplets(coef.begin(), coef.end()); + //////////////////////////////////////////////////////////////////////////// +} + diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..00e782f 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,15 +1,86 @@ #include "fd_interpolate.h" +#include -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 & W) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// +typedef Eigen::Triplet T; + +void intervalue(std::vector &coef, double aidx, double bidx, double value) { + + for (int i = 0; i < 3; i++) { + T tmp(aidx * 3 + i, bidx * 3 + i, value); + coef.push_back(tmp); + } +} + +double gidx(double i, double j, double k, int nx, int ny, int nz) { + double idx = i + nx * (j + k * ny); + assert(idx < (nx - 1) * (ny - 1) * (nz - 1)); + return idx; } + +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 & W) { + //////////////////////////////////////////////////////////////////////////// + // Add your code here + + int pnum = P.rows(); + + // size of w should be pnum * 3 x (nx-1)(ny-1)(nz-1) * 3 + std::vector coef; + for (int p = 0; p < pnum; p++) { + double px = P(p, 0) - corner[0]; + double py = P(p, 1) - corner[1]; + double pz = P(p, 2) - corner[2]; + + double xind = px / h; + double yind = py / h; + double zind = pz / h; + + // 8 equations + double xfloor = std::floor(px); + double yfloor = std::floor(py); + double zfloor = std::floor(pz); + double xceil = xfloor + 1; + double yceil = yfloor + 1; + double zceil = zfloor + 1; + + // left, left, left + double idx = gidx(xfloor, yfloor, zfloor, nx, ny, nz); + double value = (xceil - px) * (yceil - py) * (zceil - pz); + intervalue(coef, p, idx, value); + + idx = gidx(xfloor, yceil, zfloor, nx, ny, nz); + value = (xceil - px) * (py - yfloor) * (zceil - pz); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yceil, zfloor, nx, ny, nz); + value = (px - xfloor) * (py - yfloor) * (zceil - pz); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yfloor, zfloor, nx, ny, nz); + value = (px - xfloor) * (yceil - py) * (zceil - pz); + intervalue(coef, p, idx, value); + + /////////////////////////////////////////////////////// + idx = gidx(xfloor, yfloor, zceil, nx, ny, nz); + value = (xceil - px) * (yceil - px) * (pz - zfloor); + intervalue(coef, p, idx, value); + + idx = gidx(xfloor, yceil, zceil, nx, ny, nz); + value = (xceil - px) * (px - yfloor) * (pz - zfloor); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yceil, zceil, nx, ny, nz); + value = (px - xfloor) * (px - yfloor) * (pz - zfloor); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yfloor, zceil, nx, ny, nz); + value = (px - xfloor) * (yceil - py) * (pz - zfloor); + intervalue(coef, p, idx, value); + } + + W.setFromTriplets(coef.begin(), coef.end()); + + //////////////////////////////////////////////////////////////////////////// +} + diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..7762dab 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,14 +1,38 @@ #include "fd_partial_derivative.h" -void fd_partial_derivative( - const int nx, - const int ny, - const int nz, - const double h, - const int dir, - Eigen::SparseMatrix & D) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// +typedef Eigen::Triplet T; +extern int gidx(double i, double j, double k, int nx, int ny, int nz); +extern void intervalue2(std::vector &coef, int row, int aidx, int bidx); + +void fd_partial_derivative(const int nx, const int ny, const int nz, + const double h, const int dir, Eigen::SparseMatrix & D) { + //////////////////////////////////////////////////////////////////////////// + // Add your code here + // first order + std::vector coef; + for (int i = 0; i < nx - 1; i++) + for (int j = 0; j < ny - 1; j++) + for (int k = 0; k < nz - 1; k++) { + + int st = gidx(i, j, k, nx, ny, nz); + int edx = gidx(i + 1, j, k, nx, ny, nz); + int edy = gidx(i, j + 1, k, nx, ny, nz); + int edz = gidx(i, j, k + 1, nx, ny, nz); + + int row = st; + switch (dir) { + case 0: + intervalue2(coef, row, st, edx); + break; + case 1: + intervalue2(coef, row + 1, st, edy); + break; + case 2: + intervalue2(coef, row + 2, st, edz); + break; + } + } + + D.setFromTriplets(coef.begin(), coef.end()); + //////////////////////////////////////////////////////////////////////////// } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..222d2b3 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -2,55 +2,89 @@ #include #include -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); - - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// - - //////////////////////////////////////////////////////////////////////////// - // 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); +#include "fd_interpolate.h" +#include "fd_partial_derivative.h" +#include "fd_grad.h" + +#include + +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); + + //////////////////////////////////////////////////////////////////////////// + // Add your code here + //////////////////////////////////////////////////////////////////////////// + // first, we calculate spaese mtx M, which M_p3_xyz3 * X_(x-1)(y-1)(z-1)3_1 = P_p3_1 + int pnum = n; + int gnum = (nx - 1) * (ny - 1) * (nz - 1); + Eigen::SparseMatrix W(pnum * 3, gnum * 3); + fd_interpolate(nx, ny, nz, h, corner, P, W); + + // next, we calculate sparse mtx G, which G_(x-1)(y-1)(z-1)3_xyz * X_xyz_1 = D_(x-1)(y-1)(z-1)_1 + Eigen::SparseMatrix G((nx - 1) * (ny - 1) * (nz - 1) * 3, + nx * ny * nz * 1); + fd_grad(nx, ny, nz, h, G); + + // next, we combine G with M + Eigen::SparseMatrix A = W * G; + + // last, we solve Ax = B + Eigen::BiCGSTAB > solver; + solver.compute(A); + + Eigen::VectorXd b(n * 3); + for (int i = 0; i < pnum; i++) { + for (int j = 0; j < 3; j++) { + b(i * 3 + j, 0) = N(i, j); + } + } + g = solver.solve(b); + + double sigma = (W * g).mean(); + for (int i = 0; i < g.cols(); i++) + g(i) = g(i) - 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); }