From 34c4b6861746e0e6058a20c45d000bd73cdffbc9 Mon Sep 17 00:00:00 2001 From: eriszhang Date: Thu, 27 Sep 2018 11:48:53 -0400 Subject: [PATCH] finish HW1... --- src/fd_grad.cpp | 14 ++++-- src/fd_interpolate.cpp | 46 +++++++++++++++++-- src/fd_partial_derivative.cpp | 45 +++++++++++++++++-- src/poisson_surface_reconstruction.cpp | 62 +++++++++++++++++++------- 4 files changed, 142 insertions(+), 25 deletions(-) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..07cb3f0 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,4 +1,9 @@ #include "fd_grad.h" +#include "igl/cat.h" +#include "fd_partial_derivative.h" +#include + +using namespace Eigen; void fd_grad( const int nx, @@ -7,7 +12,10 @@ void fd_grad( const double h, Eigen::SparseMatrix & G) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + SparseMatrix Dx, Dy, Dz, buffer; + 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, buffer); + igl::cat(1, buffer, Dz, G); } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..80d9942 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,4 +1,8 @@ #include "fd_interpolate.h" +#include + +using namespace std; +using namespace Eigen; void fd_interpolate( const int nx, @@ -9,7 +13,43 @@ void fd_interpolate( const Eigen::MatrixXd & P, Eigen::SparseMatrix & W) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + vector> triplets; + triplets.reserve(P.rows() * 8); + + for (int i = 0; i < P.rows(); i++) { + RowVector3d diff = P.row(i) - corner; + + double xr = diff(0) * 1.0 / h; + double yr = diff(1) * 1.0 / h; + double zr = diff(2) * 1.0 / h; + + int xi = floor(xr); + double xd = xr - xi; + + int yi = floor(yr); + double yd = yr - yi; + + int zi = floor(zr); + double zd = zr - zi; + + // xi, yi, zi + triplets.push_back(Triplet(i, nx * ny * zi + nx * yi + xi, (1 - xd) * (1 - yd) * (1 - zd))); + // xi + 1, yi, zi + triplets.push_back(Triplet(i, nx * ny * zi + nx * yi + (xi + 1), xd * (1 - yd) * (1 - zd))); + // xi, yi + 1, zi + triplets.push_back(Triplet(i, nx * ny * zi + nx * (yi + 1) + xi, (1 - xd) * yd * (1 - zd))); + // xi, yi, zi + 1 + triplets.push_back(Triplet(i, nx * ny * (zi + 1) + nx * yi + xi, (1 - xd) * (1 - yd) * zd)); + // xi + 1, yi + 1, zi + triplets.push_back(Triplet(i, nx * ny * zi + nx * (yi + 1) + (xi + 1), xd * yd * (1 - zd))); + // xi + 1, yi, zi + 1 + triplets.push_back(Triplet(i, nx * ny * (zi + 1) + nx * yi + (xi + 1), xd * (1 - yd) * zd)); + // xi, yi + 1, zi + 1 + triplets.push_back(Triplet(i, nx * ny * (zi + 1) + nx * (yi + 1) + xi, (1 - xd) * yd * zd)); + // xi + 1, yi + 1, zi + 1 + triplets.push_back(Triplet(i, nx * ny * (zi + 1) + nx * (yi + 1) + (xi + 1), xd * yd * zd)); + } + + W.resize(P.rows(), nx*ny*nz); + W.setFromTriplets(triplets.begin(), triplets.end()); } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..e7d193b 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,5 +1,8 @@ #include "fd_partial_derivative.h" +using namespace Eigen; +using namespace std; + void fd_partial_derivative( const int nx, const int ny, @@ -8,7 +11,43 @@ void fd_partial_derivative( const int dir, Eigen::SparseMatrix & D) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + int nxd = nx; + int nyd = ny; + int nzd = nz; + + if (dir == 0) { + nxd--; + } + else if (dir == 1) { + nyd--; + } + else if (dir == 2) { + nzd--; + } + + vector> triplets; + triplets.reserve(D.rows()*2); + + for (int i = 0; i < nxd; i++) { + for (int j = 0; j < nyd; j++) { + for (int k = 0; k < nzd; k++) { + // same for all the three cases + int cur = nxd * nyd * k + nxd * j + i; + triplets.push_back(Triplet(cur, nx * ny * k + nx * j + i, - 1.0)); + + if (dir == 0) { + triplets.push_back(Triplet(cur, nx * ny * k + nx * j + (i + 1), 1.0)); + } + else if (dir == 1) { + triplets.push_back(Triplet(cur, nx * ny * k + nx * (j + 1) + i, 1.0)); + } + else if (dir == 2) { + triplets.push_back(Triplet(cur, nx * ny * (k + 1) + nx * j + i, 1.0)); + } + } + } + } + + D.resize(nxd*nyd*nzd, nx*ny*nz); + D.setFromTriplets(triplets.begin(), triplets.end()); } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..0b3d4a5 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,12 +1,18 @@ #include "poisson_surface_reconstruction.h" #include #include +#include "fd_grad.h" +#include "fd_interpolate.h" +#include + +using namespace std; +using namespace Eigen; void poisson_surface_reconstruction( - const Eigen::MatrixXd & P, - const Eigen::MatrixXd & N, - Eigen::MatrixXd & V, - Eigen::MatrixXi & F) + const MatrixXd & P, + const MatrixXd & N, + MatrixXd & V, + MatrixXi & F) { //////////////////////////////////////////////////////////////////////////// // Construct FD grid, CONGRATULATIONS! You get this for free! @@ -23,14 +29,14 @@ void poisson_surface_reconstruction( // 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.); + RowVector3d corner = P.colwise().minCoeff().array()-pad*h; + // Grid dimensions should be at least 3 + nx = max((P.col(0).maxCoeff()-P.col(0).minCoeff()+(2.*pad)*h)/h,3.); + ny = max((P.col(1).maxCoeff()-P.col(1).minCoeff()+(2.*pad)*h)/h,3.); + nz = 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++) + MatrixXd x(nx*ny*nz, 3); + for(int i = 0; i < nx; i++) { for(int j = 0; j < ny; j++) { @@ -38,15 +44,39 @@ void poisson_surface_reconstruction( { // Convert subscript to index const auto ind = i + nx*(j + k * ny); - x.row(ind) = corner + h*Eigen::RowVector3d(i,j,k); + x.row(ind) = corner + h * RowVector3d(i,j,k); } } } - Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz); + VectorXd g = VectorXd::Zero(nx*ny*nz); + + VectorXd vx, vy, vz, v; + SparseMatrix Wx, Wy, Wz; + //notice the shift of corner here + fd_interpolate(nx - 1, ny, nz, h, corner + RowVector3d(0.5 * h, 0, 0), P, Wx); + vx = Wx.transpose() * N.col(0); + fd_interpolate(nx, ny - 1, nz, h, corner + RowVector3d(0, 0.5 * h, 0), P, Wy); + vy = Wy.transpose() * N.col(1); + fd_interpolate(nx, ny, nz - 1, h, corner + RowVector3d(0, 0, 0.5 * h), P, Wz); + vz = Wz.transpose() * N.col(2); + + v.resize(vx.rows() + vy.rows() + vz.rows()); // remember to initialize size of v here + v << vx, vy, vz; + + SparseMatrix G, left, W; // product of sparsematrices are also sparse + VectorXd right; + fd_grad(nx, ny, nz, h, G); + left = G.transpose() * G; + right = G.transpose() * v; + BiCGSTAB> solver; + solver.compute(left); + g = solver.solve(right); + + fd_interpolate(nx, ny, nz, h, corner, P, W); + + // compute sigma and then shift g + g = (g.array() - (MatrixXd::Ones(1, n) * W * g).value() * 1.0 / n).matrix(); - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// // Run black box algorithm to compute mesh from implicit function: this