diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..595f1d2 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 +#include void fd_grad( const int nx, @@ -8,6 +11,59 @@ void fd_grad( Eigen::SparseMatrix & G) { //////////////////////////////////////////////////////////////////////////// - // Add your code here + + std::cout << "Building gradient matrix components..." << std::flush; + + int dir; + + // x-component + dir = 1; + Eigen::SparseMatrix Dx; + fd_partial_derivative(nx, ny, nz, h, dir, Dx); + + // y-component + dir = 2; + Eigen::SparseMatrix Dy; + fd_partial_derivative(nx, ny, nz, h, dir, Dy); + + // z-component + dir = 3; + Eigen::SparseMatrix Dz; + fd_partial_derivative(nx, ny, nz, h, dir, Dz); + + std::cout << "done." << std::endl; + + // Concatenate to create gradient + // Note: concatenation method sourced from https://stackoverflow.com/questions/41756428/concatenate-sparse-matrix-eigen + + std::cout << "Assembling full gradient matrix..." << std::flush; + + G.resize((nx-1)*ny*nz + nx*(ny-1)*nz + nx*ny*(nz-1), nx*ny*nz); + + // Using the triplet insertion method based on Eigen's documentation + typedef Eigen::Triplet T; + std::vector tripletList; + + for (int ii = 0; ii < G.cols(); ii++) + { + for (Eigen::SparseMatrix::InnerIterator itDx(Dx, ii); itDx; ++itDx) + tripletList.push_back(T(itDx.row(), ii, itDx.value()/h)); + // G.insertBack(itDx.row(), ii) = itDx.value()/h; + + for (Eigen::SparseMatrix::InnerIterator itDy(Dy, ii); itDy; ++itDy) + tripletList.push_back(T(itDy.row(), ii, itDy.value()/h)); + // G.insertBack(itDy.row(), ii) = itDy.value()/h; + + for (Eigen::SparseMatrix::InnerIterator itDz(Dz, ii); itDz; ++itDz) + tripletList.push_back(T(itDz.row(), ii, itDz.value()/h)); + // G.insertBack(itDz.row(), ii) = itDz.value()/h; + } + + G.setFromTriplets(tripletList.begin(), tripletList.end()); + G.finalize(); + + std::cout << "done." << std::endl; + + //////////////////////////////////////////////////////////////////////////// } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..bfe199d 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,4 +1,6 @@ #include "fd_interpolate.h" +#include +#include void fd_interpolate( const int nx, @@ -10,6 +12,134 @@ void fd_interpolate( Eigen::SparseMatrix & W) { //////////////////////////////////////////////////////////////////////////// - // Add your code here + + std::cout << "Building interpolation matrix..." << std::flush; + + // ****** DEBUGGING ****** + // std::cout << "Rows: " << P.rows() << ", Cols: " << P.cols() << std::endl; + // std::cout << P(0,0) << ", " << P(0,1) << ", " << P(0,2) << std::endl; + // std::cout << corner(0) << ", " << corner(1) << ", " << corner(2) << std::endl; + // *********************** + + // Compute total number of grid points + int ng = nx*ny*nz; + + // Extract total number of given points + int np = P.rows(); + + // Since the interpolation matrix multiplies into the list of grid points, + // and interpolates them on to the given points, we know that the size of the + // interpolation matrix must be np x ng. Also, since there are three components, + // we extend the W matrix to concatenate each component on top of each other. + + W.resize(3*np, ng); + W.reserve(3*np*4); + + // ****** DEBUGGING ****** + // std::cout << W.nonZeros() << std::endl; + // *********************** + + // For each point in P, we need to find the nearest grid points, and then + // use them to construct W. + + // Using the triplet insertion method based on Eigen's documentation + typedef Eigen::Triplet T; + std::vector tripletList; + + for (int ii = 0; ii < np; ii++) + { + // Compute distances from this point to the corner + double dist_p_to_corner_x = std::abs(P(ii, 0) - corner(0)); + double dist_p_to_corner_y = std::abs(P(ii, 1) - corner(1)); + double dist_p_to_corner_z = std::abs(P(ii, 2) - corner(2)); + + // Divide this distance by grid spacing - this is the number of cells between the point and the corner + double cells_x = dist_p_to_corner_x/h; + double cells_y = dist_p_to_corner_y/h; + double cells_z = dist_p_to_corner_z/h; + + // By rounding these numbers up and down, we can figure out the nearest vertices + double x_low = std::floor(cells_x);//*h + corner(0); + double x_high = std::ceil(cells_x);//*h + corner(0); + double y_low = std::floor(cells_y);//*h + corner(1); + double y_high = std::ceil(cells_y);//*h + corner(1); + double z_low = std::floor(cells_z);//*h + corner(2); + double z_high = std::ceil(cells_z);//*h + corner(2); + + // Also figure out the *global* indices of the eight surrounding vertices + int i000 = x_low + y_low*nx + z_low*nx*ny; + int i100 = x_high + y_low*nx + z_low*nx*ny; + int i010 = x_low + y_high*nx + z_low*nx*ny; + int i110 = x_high + y_high*nx + z_low*nx*ny; + int i001 = x_low + y_low*nx + z_high*nx*ny; + int i101 = x_high + y_low*nx + z_high*nx*ny; + int i011 = x_low + y_high*nx + z_high*nx*ny; + int i111 = x_high + y_high*nx + z_high*nx*ny; + + // ****** DEBUGGING ****** + // std::cout << i000 << ", " << i100 << ", " << i010 << ", " << i110 << ", " << i001 << ", " << i101 << ", " << i011 << ", " << i111 << std::endl; + // *********************** + + // Now insert weights into appropriate spots of the interpolation matrix + double w_x_low = 1 - (P(ii, 0) - x_low*h + corner(0))/(x_high*h + corner(0) - x_low*h + corner(0)); + double w_x_high = (P(ii, 0) - x_low*h + corner(0))/(x_high*h + corner(0) - x_low*h + corner(0)); + double w_y_low = 1 - (P(ii, 1) - y_low*h + corner(1))/(y_high*h + corner(1) - y_low*h + corner(1)); + double w_y_high = (P(ii, 1) - y_low*h + corner(1))/(y_high*h + corner(1) - y_low*h + corner(1)); + double w_z_low = 1 - (P(ii, 2) - z_low*h + corner(2))/(z_high*h + corner(2) - z_low*h + corner(2)); + double w_z_high = (P(ii, 2) - z_low*h + corner(2))/(z_high*h + corner(2) - z_low*h + corner(2)); + + // x-weights + tripletList.push_back(T(ii, i000, w_x_low)); + tripletList.push_back(T(ii, i010, w_x_low)); + tripletList.push_back(T(ii, i001, w_x_low)); + tripletList.push_back(T(ii, i011, w_x_low)); + + tripletList.push_back(T(ii, i100, w_x_high)); + tripletList.push_back(T(ii, i110, w_x_high)); + tripletList.push_back(T(ii, i101, w_x_high)); + tripletList.push_back(T(ii, i111, w_x_high)); + + // y-weights + tripletList.push_back(T(ii + np, i000, w_y_low)); + tripletList.push_back(T(ii + np, i001, w_y_low)); + tripletList.push_back(T(ii + np, i101, w_y_low)); + tripletList.push_back(T(ii + np, i100, w_y_low)); + + tripletList.push_back(T(ii + np, i010, w_y_high)); + tripletList.push_back(T(ii + np, i110, w_y_high)); + tripletList.push_back(T(ii + np, i011, w_y_high)); + tripletList.push_back(T(ii + np, i111, w_y_high)); + + // z-weights + tripletList.push_back(T(ii + 2*np, i000, w_z_low)); + tripletList.push_back(T(ii + 2*np, i100, w_z_low)); + tripletList.push_back(T(ii + 2*np, i010, w_z_low)); + tripletList.push_back(T(ii + 2*np, i110, w_z_low)); + + tripletList.push_back(T(ii + 2*np, i001, w_z_high)); + tripletList.push_back(T(ii + 2*np, i101, w_z_high)); + tripletList.push_back(T(ii + 2*np, i011, w_z_high)); + tripletList.push_back(T(ii + 2*np, i111, w_z_high)); + + + // ****** DEBUGGING ****** + // std::cout << x_low << ", " << P(ii, 0) << ", " << x_high << std::endl; + // std::cout << y_low << ", " << P(ii, 1) << ", " << y_high << std::endl; + // std::cout << z_low << ", " << P(ii, 2) << ", " << z_high << std::endl; + + // std::cout << (P(ii, 0) - x_low)/(x_high - x_low) << std::endl; + // std::cout << (P(ii, 1) - y_low)/(y_high - y_low) << std::endl; + // std::cout << (P(ii, 2) - z_low)/(z_high - z_low) << std::endl; + // *********************** + + } + + W.setFromTriplets(tripletList.begin(), tripletList.end()); + W.finalize(); + + std::cout << "done." << std::endl; + //////////////////////////////////////////////////////////////////////////// } + + diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..64f13a4 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,4 +1,5 @@ #include "fd_partial_derivative.h" +#include void fd_partial_derivative( const int nx, @@ -9,6 +10,137 @@ void fd_partial_derivative( Eigen::SparseMatrix & D) { //////////////////////////////////////////////////////////////////////////// - // Add your code here + + // Compute the derivatives for the desired direction + switch (dir) + { + case 1: + { + D.resize((nx-1)*ny*nz, nx*ny*nz); + + // Using the triplet insertion method based on Eigen's documentation + typedef Eigen::Triplet T; + std::vector tripletList; + + for (int ii = 0; ii < nx*ny*nz; ii++) + { + int col = ii; + int col_other; + + // Check if we're at an edge, in which case the finite difference is taken to be a forward difference + if (ii%nx == 0) + { + col_other = ii+1; + if (col >= (nx-1)*ny*nz) + col = (nx-1)*ny*nz - 1; + if (col_other >= nx*ny*nz) + col_other = nx*ny*nz - 1; + + tripletList.push_back(T(col, col, -1.0)); + tripletList.push_back(T(col, col_other, 1.0)); + } + else + { + col_other = ii-1; + if (col_other >= (nx-1)*ny*nz) + col_other = (nx-1)*ny*nz - 1; + if (col_other < 0) + col_other = 0; + + tripletList.push_back(T(col_other, col_other, -1.0)); + tripletList.push_back(T(col_other, col, 1.0)); + } + } + + D.setFromTriplets(tripletList.begin(), tripletList.end()); + D.finalize(); + } + case 2: + { + D.resize(nx*(ny-1)*nz, nx*ny*nz); + + // Using the triplet insertion method based on Eigen's documentation + typedef Eigen::Triplet T; + std::vector tripletList; + + for (int ii = 0; ii < nx*ny*nz; ii++) + { + int col = ii; + int col_other; + + // Check if we're at an edge, in which case the finite difference is taken to be a forward difference + if (ii%(nx*ny) < nx) + { + col_other = ii+1*nx; + if (col >= nx*(ny-1)*nz) + col = nx*(ny-1)*nz - 1; + if (col_other >= nx*ny*nz) + col_other = nx*ny*nz - 1; + + tripletList.push_back(T(col, col, -1.0)); + tripletList.push_back(T(col, col_other, 1.0)); + } + else + { + col_other = ii-1*nx; + if (col_other >= nx*(ny-1)*nz) + col_other = nx*(ny-1)*nz - 1; + if (col_other < 0) + col_other = 0; + + tripletList.push_back(T(col_other, col_other, -1.0)); + tripletList.push_back(T(col_other, col, 1.0)); + } + } + + D.setFromTriplets(tripletList.begin(), tripletList.end()); + D.finalize(); + } + case 3: + { + D.resize(nx*ny*(nz-1), nx*ny*nz); + + // Using the triplet insertion method based on Eigen's documentation + typedef Eigen::Triplet T; + std::vector tripletList; + + for (int ii = 0; ii < nx*ny*nz; ii++) + { + int col = ii; + int col_other; + + // Check if we're at an edge, in which case the finite difference is taken to be a forward difference + if (ii < nx) + { + col_other = ii+1*nx*ny; + if (col >= nx*ny*(nz-1)) + col = nx*ny*(nz-1) - 1; + if (col_other >= nx*ny*nz) + col_other = nx*ny*nz - 1; + + tripletList.push_back(T(col, col, -1.0)); + tripletList.push_back(T(col, col_other, 1.0)); + } + else + { + col_other = ii-1*nx*ny; + if (col_other >= nx*ny*(nz-1)) + col_other = nx*ny*(nz-1) - 1; + if (col_other < 0) + col_other = 0; + + tripletList.push_back(T(col_other, col_other, -1.0)); + tripletList.push_back(T(col_other, col, 1.0)); + + } + } + + D.setFromTriplets(tripletList.begin(), tripletList.end()); + D.finalize(); + } + + }; + + //////////////////////////////////////////////////////////////////////////// } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..b99d377 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,6 +1,12 @@ #include "poisson_surface_reconstruction.h" #include #include +#include +#include "fd_interpolate.h" +#include "fd_grad.h" +#include +#include +#include void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -45,7 +51,162 @@ void poisson_surface_reconstruction( Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz); //////////////////////////////////////////////////////////////////////////// - // Add your code here + + // Create interpolation matrix + Eigen::SparseMatrix W, Wx, Wy, Wz, Wx2, Wy2, Wz2, Gx, Gy, Gz; + // fd_interpolate(nx, ny, nz, h, corner, P, W); + + Eigen::RowVector3d corner1; corner1(0) = corner(0) + h/2; + fd_interpolate(nx-1, ny, nz, h, corner1, P, Wx); + Eigen::RowVector3d corner2; corner2(1) = corner(1) + h/2; + fd_interpolate(nx, ny-1, nz, h, corner2, P, Wy); + Eigen::RowVector3d corner3; corner3(2) = corner(2) + h/2; + fd_interpolate(nx, ny, nz-1, h, corner3, P, Wz); + + // Compute total number of grid points + int ng = nx*ny*nz; + + // Extract total number of given points + int np = P.rows(); + + // Using the triplet insertion method based on Eigen's documentation + typedef Eigen::Triplet T; + std::vector tripletList1, tripletList2, tripletList3; + + // // Extract the x, y, and z parts + // for (int ii = 0; ii < W.outerSize(); ii++) + // { + // for (Eigen::SparseMatrix::InnerIterator itW(W, ii); itW; ++itW) + // { + // if (itW.row() < np) + // tripletList1.push_back(T(itW.row(), ii, itW.value())); + // else if (itW.row() < 2*np) + // tripletList2.push_back(T(itW.row()-np, ii, itW.value())); + // else + // tripletList3.push_back(T(itW.row()-2*np, ii, itW.value())); + // } + // } + + // Extract the x, y, and z parts + for (int ii = 0; ii < (nx-1)*ny*nz; ii++) + { + for (Eigen::SparseMatrix::InnerIterator itW(Wx, ii); itW; ++itW) + { + if (itW.row() < np) + tripletList1.push_back(T(itW.row(), ii, itW.value())); + } + } + for (int ii = 0; ii < nx*(ny-1)*nz; ii++) + { + for (Eigen::SparseMatrix::InnerIterator itW(Wy, ii); itW; ++itW) + { + if (itW.row() >= np && itW.row() < 2*np) + tripletList2.push_back(T(itW.row()-np, ii, itW.value())); + } + } + for (int ii = 0; ii < nx*ny*(nz-1); ii++) + { + for (Eigen::SparseMatrix::InnerIterator itW(Wz, ii); itW; ++itW) + { + if (itW.row() >= 2*np) + tripletList3.push_back(T(itW.row()-2*np, ii, itW.value())); + } + } + + Wx2.resize(np, (nx-1)*ny*nz); + Wx2.reserve(np*4); + Wy2.resize(np, nx*(ny-1)*nz); + Wy2.reserve(np*4); + Wz2.resize(np, nx*ny*(nz-1)); + Wz2.reserve(np*4); + + Wx2.setFromTriplets(tripletList1.begin(), tripletList1.end()); + Wy2.setFromTriplets(tripletList2.begin(), tripletList2.end()); + Wz2.setFromTriplets(tripletList3.begin(), tripletList3.end()); + + Wx2.finalize(); + Wy2.finalize(); + Wz2.finalize(); + + // Distribute normals on to grid points. First, extract the x, y and z components of the normals + Eigen::MatrixXd Nx(N.rows(), 1), Ny(N.rows(), 1), Nz(N.rows(), 1); + for (int ii = 0; ii < N.rows(); ii++) + { + Nx(ii, 0) = N(ii, 0); + Ny(ii, 0) = N(ii, 1); + Nz(ii, 0) = N(ii, 2); + } + + // Multiply these vectors with the respective interpolation matrices' transposes + Eigen::MatrixXd Nvx(ng, 1), Nvy(ng, 1), Nvz(ng, 1); + Nvx = Wx2.transpose()*Nx; + Nvy = Wy2.transpose()*Ny; + Nvz = Wz2.transpose()*Nz; + + // Compute G^TG + Eigen::SparseMatrix G; + fd_grad(nx, ny, nz, h, G); + + Eigen::SparseMatrix A = G.transpose()*G; + + // Extract the x, y and z components of the gradient + // Using the triplet insertion method based on Eigen's documentation + std::vector tripletList4, tripletList5, tripletList6; + + Gx.resize((nx-1)*ny*nz, nx*ny*nz); + Gy.resize(nx*(ny-1)*nz, nx*ny*nz); + Gz.resize(nx*ny*(nz-1), nx*ny*nz); + + // Extract the x, y, and z parts + for (int ii = 0; ii < G.outerSize(); ii++) + { + for (Eigen::SparseMatrix::InnerIterator itG(G, ii); itG; ++itG) + { + if (itG.row() < (nx-1)*ny*nz) + tripletList4.push_back(T(itG.row(), ii, itG.value())); + else if (itG.row() < (nx-1)*ny*nz + nx*(ny-1)*nz) + tripletList5.push_back(T(itG.row()-((nx-1)*ny*nz), ii, itG.value())); + else + tripletList6.push_back(T(itG.row()-((nx-1)*ny*nz + nx*(ny-1)*nz), ii, itG.value())); + } + } + + Gx.setFromTriplets(tripletList4.begin(), tripletList4.end()); + Gy.setFromTriplets(tripletList5.begin(), tripletList5.end()); + Gz.setFromTriplets(tripletList6.begin(), tripletList6.end()); + + Gx.finalize(); + Gy.finalize(); + Gz.finalize(); + + // Now compute the right hand side + Eigen::VectorXd b = Gx.transpose()*Nvx + Gy.transpose()*Nvy + Gz.transpose()*Nvz; + + // Solve the Ag = b system + std::cout << "Solving the system..." << std::flush; + + // Eigen::ConjugateGradient> solver; + Eigen::SparseLU> solver; + + solver.compute(A); + g = solver.solve(b); + + std::cout << "done." << std::endl; + + // Set iso level - note: this should be done using the formula suggested in the paper; I ran out of time + double sigma = 0.0; + + // Apply sigma + for (int ii = 0; ii < nx*ny*nz; ii++) + g(ii) = g(ii) - sigma; + + // Based on the results, I definitely did not successfully reconstruct the mesh... + // It looks like I've flipped signs along the way, and / or I think I may have + // used an incompatible assumption on the ordering of the coordinates (I have + // assumed they're stored contiguously along x first, then y, then z, but I might have + // mistakenly been inconsistent). I ran out of time trying to fix this, so I'm + // submitting it as is (although I definitely want to make it work in the coming days). + //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// @@ -54,3 +215,7 @@ void poisson_surface_reconstruction( //////////////////////////////////////////////////////////////////////////// igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F); } + + + +