From 28d97502a0c284ed5ed41af6afe59057501ed276 Mon Sep 17 00:00:00 2001 From: zero-impact Date: Mon, 29 Jan 2018 20:56:15 -0500 Subject: [PATCH] Done. --- include/poisson_surface_reconstruction.h | 2 +- src/fd_grad.cpp | 18 ++++++++-- src/fd_interpolate.cpp | 46 ++++++++++++++++++++++-- src/fd_partial_derivative.cpp | 42 ++++++++++++++++++++-- src/poisson_surface_reconstruction.cpp | 37 +++++++++++++++++++ 5 files changed, 135 insertions(+), 10 deletions(-) diff --git a/include/poisson_surface_reconstruction.h b/include/poisson_surface_reconstruction.h index c0ce15d..fac6fe1 100644 --- a/include/poisson_surface_reconstruction.h +++ b/include/poisson_surface_reconstruction.h @@ -1,7 +1,7 @@ #ifndef POISSON_SURFACE_RECONSTRUCTION_H #define POISSON_SURFACE_RECONSTRUCTION_H #include - +#include // Takes input sample points P and input normals N // and gives a watertight mesh using a simplified // version of [Kazhdan et. al 2006] diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..a85c1cf 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,4 +1,7 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" + +#include "igl/cat.h" void fd_grad( const int nx, @@ -7,7 +10,16 @@ void fd_grad( const double h, Eigen::SparseMatrix & G) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + Eigen::SparseMatrix Dx; + fd_partial_derivative(nx, ny, nz, h, 0, Dx); + + Eigen::SparseMatrix Dy; + fd_partial_derivative(nx, ny, nz, h, 1, Dy); + + Eigen::SparseMatrix Dz; + fd_partial_derivative(nx, ny, nz, h, 2, Dz); + + Eigen::SparseMatrix T; + igl::cat(1, Dx, Dy, T); + igl::cat(1, T, Dz, G); } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..5ebdc4e 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,5 +1,10 @@ #include "fd_interpolate.h" +#include +#include + +#include "igl/floor.h" + void fd_interpolate( const int nx, const int ny, @@ -9,7 +14,42 @@ void fd_interpolate( const Eigen::MatrixXd & P, Eigen::SparseMatrix & W) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + // Get the grid corner indices for each input point + Eigen::MatrixXd P_rel_grid = (P.rowwise() - corner) / h; + + Eigen::MatrixXi P_corners_ijk; + igl::floor(P_rel_grid, P_corners_ijk); + + // These are the first three weights + Eigen::MatrixXd P_rel_corner = P_rel_grid - P_corners_ijk.cast(); + + auto ind = [&nx, &ny](int i, int j, int k) { + return i + nx * (j + k * ny); + }; + + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve(8 * P.rows()); + for(int l = 0; l < P.rows(); l++) + { + int i = P_corners_ijk(l, 0); + int j = P_corners_ijk(l, 1); + int k = P_corners_ijk(l, 2); + + double wx = P_rel_corner(l, 0); + double wy = P_rel_corner(l, 1); + double wz = P_rel_corner(l, 2); + + tripletList.push_back(T(l, ind(i+0, j+0, k+0), (1. - wx) * (1. - wy) * (1. - wz))); + tripletList.push_back(T(l, ind(i+0, j+0, k+1), (1. - wx) * (1. - wy) * wz)); + tripletList.push_back(T(l, ind(i+0, j+1, k+0), (1. - wx) * wy * (1. - wz))); + tripletList.push_back(T(l, ind(i+0, j+1, k+1), (1. - wx) * wy * wz)); + tripletList.push_back(T(l, ind(i+1, j+0, k+0), wx * (1. - wy) * (1. - wz))); + tripletList.push_back(T(l, ind(i+1, j+0, k+1), wx * (1. - wy) * wz)); + tripletList.push_back(T(l, ind(i+1, j+1, k+0), wx * wy * (1. - wz))); + tripletList.push_back(T(l, ind(i+1, j+1, k+1), wx * wy * wz)); + } + + W.resize(P.rows(), nx*ny*nz); + W.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..1437828 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -8,7 +8,43 @@ void fd_partial_derivative( const int dir, Eigen::SparseMatrix & D) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + int is_x = dir == 0; + int is_y = dir == 1; + int is_z = dir == 2; + + int nxd = nx - 1 * is_x; + int nyd = ny - 1 * is_y; + int nzd = nz - 1 * is_z; + + auto ind = [&nx, &ny](int i, int j, int k) { + return i + nx * (j + k * ny); + }; + + auto ind_d = [&nxd, &nyd](int i, int j, int k) { + return i + nxd * (j + k * nyd); + }; + + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve(nxd*nyd*nzd); + + for(int i = 0; i < nxd; i++) + { + for(int j = 0; j < nyd; j++) + { + for(int k = 0; k < nzd; k++) + { + int ip = i + 1 * is_x; + int jp = j + 1 * is_y; + int kp = k + 1 * is_z; + + int row = ind_d(i, j, k); + tripletList.push_back(T(row, ind(i,j,k), 1./h)); + tripletList.push_back(T(row, ind(ip,jp,kp), -1./h)); + } + } + } + + D.resize(nxd*nyd*nzd, 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..c11d140 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,6 +1,9 @@ #include "poisson_surface_reconstruction.h" #include #include +#include +#include "fd_interpolate.h" +#include "fd_grad.h" void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -47,7 +50,41 @@ void poisson_surface_reconstruction( //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// + Eigen::RowVector3d dx(h/2., 0., 0.); + Eigen::RowVector3d dy(0., h/2., 0.); + Eigen::RowVector3d dz(0., 0., h/2.); + Eigen::SparseMatrix Wx; + fd_interpolate(nx - 1, ny, nz, h, corner + dx, P.rowwise() + dx, Wx); + Eigen::VectorXd vx = Wx.transpose() * N.col(0); + std::cout << vx.size() << std::endl; + Eigen::SparseMatrix Wy; + fd_interpolate(nx, ny - 1, nz, h, corner + dy, P.rowwise() + dy, Wy); + Eigen::VectorXd vy = Wy.transpose() * N.col(1); + + Eigen::SparseMatrix Wz; + fd_interpolate(nx, ny, nz - 1, h, corner + dz, P.rowwise() + dz, Wz); + Eigen::VectorXd vz = Wz.transpose() * N.col(2); + + Eigen::VectorXd v(vx.size() + vy.size() + vz.size()); + v << vx, vy, vz; + + Eigen::SparseMatrix G; + fd_grad(nx, ny, nz, h, G); + std::cout << "G rows: " << G.rows() << " G cols: " << G.cols() << std::endl; + + std::cout << "Solving..." << std::endl; + Eigen::ConjugateGradient> solver; + solver.compute(G.transpose() * G); + g = solver.solve(G.transpose() * v); + std::cout << "#iterations: " << solver.iterations() << std::endl; + std::cout << "estimated error: " << solver.error() << std::endl; + + Eigen::SparseMatrix W; + fd_interpolate(nx, ny, nz, h, corner, P, W); + + double sigma = (W * g).mean(); + g = -(g - Eigen::VectorXd::Ones(g.size()) * 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