diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..8c6e8e7 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,13 +1,19 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" +#include void fd_grad( - const int nx, - const int ny, - const int nz, - const double h, - Eigen::SparseMatrix & G) + const int nx, + const int ny, + const int nz, + const double h, + Eigen::SparseMatrix & G) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + Eigen::SparseMatrix Dx, Dy, Dz, Dxy; + 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); + + igl::cat(1, Dx, Dy, Dxy); + igl::cat(1, Dxy, Dz, G); } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..f74f0ae 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,15 +1,44 @@ #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) + 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; + std::vector tripletList; + tripletList.reserve(P.rows() * 8); + for(int i = 0; i < P.rows(); i++) { + // relative position to corner + double x = (P(i, 0) - corner[0])/h; + double y = (P(i, 1) - corner[1])/h; + double z = (P(i, 2) - corner[2])/h; + + // previous on grid position + int x1 = floor(x); + int y1 = floor(y); + int z1 = floor(z); + + // next on grid position + int x2 = x1 + 1; + int y2 = y1 + 1; + int z2 = z1 + 1; + + // calculate weight for each of eight surrounding point + tripletList.push_back(T(i, (x1 + y1 * nx + z1 * ny * nx), (x2 - x) * (y2 - y) * (z2 - z))); + tripletList.push_back(T(i, (x1 + y1 * nx + z2 * ny * nx), (x2 - x) * (y2 - y) * (z - z1))); + tripletList.push_back(T(i, (x1 + y2 * nx + z1 * ny * nx), (x2 - x) * (y - y1) * (z2 - z))); + tripletList.push_back(T(i, (x1 + y2 * nx + z2 * ny * nx), (x2 - x) * (y - y1) * (z - z1))); + tripletList.push_back(T(i, (x2 + y1 * nx + z1 * ny * nx), (x - x1) * (y2 - y) * (z2 - z))); + tripletList.push_back(T(i, (x2 + y1 * nx + z2 * ny * nx), (x - x1) * (y2 - y) * (z - z1))); + tripletList.push_back(T(i, (x2 + y2 * nx + z1 * ny * nx), (x - x1) * (y - y1) * (z2 - z))); + tripletList.push_back(T(i, (x2 + y2 * nx + z2 * ny * nx), (x - x1) * (y - y1) * (z - z1))); + } + W.resize(P.rows(), nx * ny * nz); + W.setFromTriplets(tripletList.begin(), tripletList.end()); +} \ No newline at end of file diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..4662260 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,14 +1,54 @@ #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) + const int nx, + const int ny, + const int nz, + const double h, + const int dir, + Eigen::SparseMatrix & D) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + // compute number of entries + int x = nx; + int y = ny; + int z = nz; + switch (dir){ + case 0: + x--; + break; + case 1: + y--; + break; + case 2: + z--; + break; + } + + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve(x * y * z * 2); + + for(int i = 0; i < x; i++) { + for (int j = 0; j < y; j++) { + for (int k = 0; k < z; k++) { + int s = i + j * x + k * y * x; + int t = i + j * nx + k * ny * nx; + tripletList.push_back(T(s, t, -1/h)); + switch (dir){ + case 0: + t = (i + 1) + j * nx + k * ny * nx; + break; + case 1: + t = i + (j + 1) * nx + k * ny * nx; + break; + case 2: + t = i + j * nx + (k + 1) * ny * nx; + break; + } + tripletList.push_back(T(s, t, 1/h)); + } + } + } + D.resize(x * y * z, nx * ny * nz); + D.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..5e610b4 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,56 +1,91 @@ #include "poisson_surface_reconstruction.h" #include #include +#include "fd_grad.h" +#include "fd_interpolate.h" +#include void poisson_surface_reconstruction( - const Eigen::MatrixXd & P, - const Eigen::MatrixXd & N, - Eigen::MatrixXd & V, - Eigen::MatrixXi & F) + 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); + //////////////////////////////////////////////////////////////////////////// + // 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 - //////////////////////////////////////////////////////////////////////////// + Eigen::SparseMatrix Wx, Wy, Wz, G, W; + // get gradient + fd_grad(nx, ny, nz, h, G); - //////////////////////////////////////////////////////////////////////////// - // 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); + // staggered grid + 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); + + // construct v + Eigen::VectorXd vx(nx*ny*nz - ny*nz), vy(nx*ny*nz - nx*nz), vz(nx*ny*nz - nx*ny), + 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); + + // solve for g + Eigen::BiCGSTAB> solver; + solver.compute(G.transpose() * G); + g = solver.solve(G.transpose() * v); + + // compute iso-value + fd_interpolate(nx, ny, nz, h, corner, P, W); + + // shift + // g -= (1.0 / n) * Eigen::RowVectorXd::Ones(n) * W * g * Eigen::VectorXd::Ones(nx * ny * nz); + double sigma = (W*g).sum()/n; + g = (g.array() - sigma).matrix(); + + //////////////////////////////////////////////////////////////////////////// + // 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); + // std::cout << "V is of size " << V.size() << std::endl; + // std::cout << "F is of size " << F.size() << std::endl; }