diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..dcf3283 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,4 +1,6 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" +#include "igl/cat.h" void fd_grad( const int nx, @@ -7,7 +9,14 @@ void fd_grad( const double h, Eigen::SparseMatrix & G) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + Eigen::SparseMatrix D1; + Eigen::SparseMatrix D2; + Eigen::SparseMatrix D3; + + fd_partial_derivative(nx, ny, nz, h, 0, D1); + fd_partial_derivative(nx, ny, nz, h, 1, D2); + fd_partial_derivative(nx, ny, nz, h, 2, D3); + + G = igl::cat(1, D1, D2); + G = igl::cat(1, G, D3); } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..7f892eb 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,4 +1,5 @@ #include "fd_interpolate.h" +#include void fd_interpolate( const int nx, @@ -9,7 +10,30 @@ void fd_interpolate( const Eigen::MatrixXd & P, Eigen::SparseMatrix & W) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + W.resize(P.rows(), nx*ny*nz); + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve(P.rows()*8); + + for(int i = 0; i < P.rows(); i += 1){ + + int x = floor((P(i, 0) - corner(0))/h); + int y = floor((P(i, 1) - corner(1))/h); + int z = floor((P(i, 2) - corner(2))/h); + + double dx = (P(i, 0) - corner(0))/h - x; + double dy = (P(i, 1) - corner(1))/h - y; + double dz = (P(i, 2) - corner(2))/h - z; + + tripletList.push_back(T(i, x + nx*y + nx*ny*z, (1-dx)*(1-dy)*(1-dz))); + tripletList.push_back(T(i, x + nx*y + nx*ny*(z+1), (1-dx)*(1-dy)*dz)); + tripletList.push_back(T(i, x + nx*(y+1) + nx*ny*z, (1-dx)*dy*(1-dz))); + tripletList.push_back(T(i, x + nx*(y+1) + nx*ny*(z+1), (1-dx)*dy*dz)); + tripletList.push_back(T(i, (x+1) + nx*y + nx*ny*z, dx*(1-dy)*(1-dz))); + tripletList.push_back(T(i, (x+1) + nx*y + nx*ny*(z+1), dx*(1-dy)*dz)); + tripletList.push_back(T(i, (x+1) + nx*(y+1) + nx*ny*z, dx*dy*(1-dz))); + tripletList.push_back(T(i, (x+1) + nx*(y+1) + nx*ny*(z+1), dx*dy*dz)); + + } + W.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..eb3283e 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,5 +1,5 @@ #include "fd_partial_derivative.h" - +#include void fd_partial_derivative( const int nx, const int ny, @@ -8,7 +8,25 @@ void fd_partial_derivative( const int dir, Eigen::SparseMatrix & D) { - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + int a = 1 - ceil(dir/2.0); + int b = ceil(dir/2.0) - floor(dir/2.0); + int c = floor(dir/2.0); + + int new_x = nx-a; + int new_y = ny-b; + int new_z = nz-c; + + int next = a + (b*nx) + (c*nx*ny); + + D.resize(new_x*new_y*new_z, nx*ny*nz); + + for(int i = 0; i < new_x; i += 1){ + for(int j = 0; j < new_y; j += 1){ + for(int k = 0; k < new_z; k += 1){ + D.insert(i + (new_x*j) + (new_x*new_y*k), i + (nx*j) + (nx*ny*k)) = -1; + D.insert(i + (new_x*j) + (new_x*new_y*k), i + (nx*j) + (nx*ny*k) + next) = 1; + } + } + } + std::cout << D.rows() << std::endl; } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..a9e2b05 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,6 +1,12 @@ #include "poisson_surface_reconstruction.h" +#include "fd_grad.h" +#include "fd_partial_derivative.h" +#include "fd_interpolate.h" #include #include +#include +#include "igl/cat.h" +#include void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -12,6 +18,7 @@ void poisson_surface_reconstruction( // Construct FD grid, CONGRATULATIONS! You get this for free! //////////////////////////////////////////////////////////////////////////// // number of input points + const int n = P.rows(); // Grid dimensions int nx, ny, nz; @@ -44,10 +51,50 @@ void poisson_surface_reconstruction( } Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz); - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + + Eigen::SparseMatrix G; + + Eigen::SparseMatrix W; + Eigen::SparseMatrix Wx; + Eigen::SparseMatrix Wy; + Eigen::SparseMatrix Wz; + + Eigen::RowVector3d shift_x(h/2, 0., 0.); + Eigen::RowVector3d shift_y(0., h/2, 0.); + Eigen::RowVector3d shift_z(0., 0., h/2); + + Eigen::VectorXd vx; + Eigen::VectorXd vy; + Eigen::VectorXd vz; + Eigen::VectorXd v; + + fd_grad(nx, ny, nz, h, G); + + fd_interpolate(nx, ny, nz, h, corner, P, W); + fd_interpolate(nx-1, ny, nz, h, corner + shift_x, P, Wx); + fd_interpolate(nx, ny-1, nz, h, corner + shift_y, P, Wy); + fd_interpolate(nx, ny, nz-1, h, corner + shift_z, P, Wz); + + vx = Wx.transpose() * N.col(0); + vy = Wy.transpose() * N.col(1); + vz = Wz.transpose() * N.col(2); + + v = igl::cat(1, vx, vy); + v = igl::cat(1, v, vz); + + + Eigen::BiCGSTAB> solver; + + solver.compute(G.transpose() * G); + std::cout << nx << std::endl; + std::cout << ny << std::endl; + std::cout << nz << std::endl; + std::cout << N.rows() << std::endl; + g = solver.solve(G.transpose() * v); + double sigma = (1/n)*Eigen::VectorXd::Ones(W.rows()).transpose()*(W*g); + g = g - (Eigen::VectorXd::Ones(g.rows())*sigma); + std::cout << sigma << std::endl; //////////////////////////////////////////////////////////////////////////// // Run black box algorithm to compute mesh from implicit function: this // function always extracts g=0, so "pre-shift" your g values by -sigma