diff --git a/include/fd_grad.h b/include/fd_grad.h index fa0df87..fa9fe22 100644 --- a/include/fd_grad.h +++ b/include/fd_grad.h @@ -1,6 +1,7 @@ #ifndef FD_GRAD_H #define FD_GRAD_H #include + // Construct a gradient matrix for a finite-difference grid // // Inputs: diff --git a/main.cpp b/main.cpp index 1845980..ea8db77 100644 --- a/main.cpp +++ b/main.cpp @@ -9,6 +9,12 @@ #include #include +#ifndef AP_WINDOWS_BUILD +#define PATH_PREFIX "../" +#else +#define PATH_PREFIX "" +#endif + int main(int argc, char *argv[]) { // Load in points + normals from .pwn file @@ -18,7 +24,7 @@ int main(int argc, char *argv[]) std::vector > vD; std::string line; std::fstream in; - in.open(argc>1?argv[1]:"../shared/data/hand.pwn"); + in.open(argc>1?argv[1]:PATH_PREFIX "shared/data/hand.pwn"); while(in) { std::getline(in, line); diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..b48c6bb 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,13 +1,50 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" + +typedef Eigen::Triplet Triplet; + +/** +* @brief Concatenates vertically. +*/ +void concatSparse(Eigen::SparseMatrix& source, + std::vector& targetTriplets, int firstTargetRow) { + typedef Eigen::SparseMatrix::InnerIterator Iterator; + + // From Eigen tutorials... + for (int k = 0; k < source.outerSize(); k++) { + for (Iterator it(source, k); it; ++it) { + int sourceRow = it.row(); + int column = it.col(); + double value = it.value(); + targetTriplets.push_back(Triplet(firstTargetRow + sourceRow, column, value)); + } + } + +} void fd_grad( - const int nx, - const int ny, - const int nz, - const double h, - Eigen::SparseMatrix & G) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + const int nx, + const int ny, + const int nz, + const double h, + Eigen::SparseMatrix & G) { + + Eigen::SparseMatrix Dx((nx - 1) * ny * nz, nx * ny * nz); + fd_partial_derivative(nx, ny, nz, h, 0, Dx); + + Eigen::SparseMatrix Dy(nx * (ny - 1) * nz, nx * ny * nz); + fd_partial_derivative(nx, ny, nz, h, 1, Dy); + + Eigen::SparseMatrix Dz(nx * ny * (nz - 1), nx * ny * nz); + fd_partial_derivative(nx, ny, nz, h, 2, Dz); + + // Concatenate + std::vector triplets; + concatSparse(Dx, triplets, 0); + concatSparse(Dy, triplets, Dx.rows()); + concatSparse(Dz, triplets, Dx.rows() + Dy.rows()); + + G.resize(Dx.rows() + Dy.rows() + Dz.rows(), Dx.cols()); + G.setFromTriplets(triplets.begin(), triplets.end()); + } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..dd89b8c 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,15 +1,61 @@ #include "fd_interpolate.h" 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) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + const int nx, + const int ny, + const int nz, + const double h, + const Eigen::RowVector3d & corner, + const Eigen::MatrixXd & P, + Eigen::SparseMatrix & W) { + W.resize(P.rows(), nx * ny * nz); + + auto toGridIndex = [nx, ny, nz](int x, int y, int z) -> int { + return x + nx * y + nx * ny * z; + }; + + // As the Eigen guide says to, we build the matrix from triplets. + typedef Eigen::Triplet Triplet; + std::vector tripletList; + + int rowCount = P.rows(); + for (int r = 0; r < rowCount; r++) { + + // How many h's from the corner the point is. + double xFromStart = (P(r, 0) - corner(0)) / h; + double yFromStart = (P(r, 1) - corner(1)) / h; + double zFromStart = (P(r, 2) - corner(2)) / h; + + // Indices of the lower corner of the grid cell. + int xGridIndex = floor(xFromStart); + int yGridIndex = floor(yFromStart); + int zGridIndex = floor(zFromStart); + + // Where within the current grid cell the point is. + double x = xFromStart - xGridIndex; + double y = yFromStart - yGridIndex; + double z = zFromStart - zGridIndex; + + // The weights computed by hand + double w0 = (1 - x) * (1 - z) * (1 - y); + double w1 = (1 - x) * z * (1 - y); + double w2 = (1 - x) * (1 - z) * y; + double w3 = (1 - x) * z * y; + double w4 = x * (1 - z) * (1 - y); + double w5 = x * z * (1 - y); + double w6 = x * (1 - z) * y; + double w7 = x * z * y; + + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex, zGridIndex), w0)); + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex, zGridIndex + 1), w1)); + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex + 1, zGridIndex), w2)); + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex + 1, zGridIndex + 1), w3)); + + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex, zGridIndex), w4)); + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex, zGridIndex + 1), w5)); + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex + 1, zGridIndex), w6)); + tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex + 1, zGridIndex + 1), w7)); + } + + W.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..fcab0b1 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,14 +1,45 @@ #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) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// + const int nx, + const int ny, + const int nz, + const double h, + const int dir, + Eigen::SparseMatrix & D) { + + // As the Eigen guide says to, we build the matrix from triplets. + typedef Eigen::Triplet Triplet; + std::vector tripletList; + + // Iterate over all cells. The limit depends on which direction we're + // going, so we can avoid going off the end. + int xLimit = dir == 0 ? nx - 1 : nx; + int yLimit = dir == 1 ? ny - 1 : ny; + int zLimit = dir == 2 ? nz - 1 : nz; + int r = 0; + for (int x = 0; x < xLimit; x++) { + for (int y = 0; y < yLimit; y++) { + for (int z = 0; z < zLimit; z++) { + + // Co-efficient corresponding to one of the points in the + // finite difference. I divide by h to represent the division + // by dx. + int firstIndex = x + y * nx + z * nx * ny; + tripletList.push_back(Triplet(r, firstIndex, -1.0 / h)); + + // Second co-efficient. Again we have to divide by h. + int secondX = dir == 0 ? x + 1 : x; + int secondY = dir == 1 ? y + 1 : y; + int secondZ = dir == 2 ? z + 1 : z; + int secondIndex = secondX + secondY * nx + + secondZ * nx * ny; + tripletList.push_back(Triplet(r, secondIndex, 1.0 / h)); + + r++; + } + } + } + + D.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..05ff9c3 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,6 +1,10 @@ #include "poisson_surface_reconstruction.h" #include #include +#include "fd_interpolate.h" +#include "fd_grad.h" +#include + void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -8,49 +12,95 @@ void poisson_surface_reconstruction( 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); - - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// - - //////////////////////////////////////////////////////////////////////////// - // 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); + //////////////////////////////////////////////////////////////////////////// + // 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); + } + } + } + + /// Implementation Note: + /// Was not able to get it to correctly form the mesh. There must be a + /// bug somewhere, but despite long hours of searching, I have not been + /// able to solve it :(. With the exception of the bug, I am very + /// confident in my implementation. + + // Weight matrices for each axis and the full weight matrix for the + // non-staggered grid. + Eigen::SparseMatrix xWeightMatrix, yWeightMatrix, zWeightMatrix, + fullWeightMatrix; + fd_interpolate(nx - 1, ny, nz, h, corner + Eigen::RowVector3d(h / 2, 0, 0), P, xWeightMatrix); + fd_interpolate(nx, ny - 1, nz, h, corner + Eigen::RowVector3d(0, h / 2, 0), P, yWeightMatrix); + fd_interpolate(nx, ny, nz - 1, h, corner + Eigen::RowVector3d(0, 0, h / 2), P, zWeightMatrix); + fd_interpolate(nx, ny, nz, h, corner, P, fullWeightMatrix); + + // Normals distributed onto grid. + Eigen::MatrixXd vX = xWeightMatrix.transpose() * N.col(0); + Eigen::MatrixXd vY = yWeightMatrix.transpose() * N.col(1); + Eigen::MatrixXd vZ = zWeightMatrix.transpose() * N.col(2); + Eigen::MatrixXd v(vX.rows() + vY.rows() + vZ.rows(), 1); + v << vX, vY, vZ; + + // Gradient matrix via finite differences. + Eigen::SparseMatrix gradient; + fd_grad(nx, ny, nz, h, gradient); + + // Solve for the grid values. The standard form for least squares is + // something like: Ax = b, so we have to compute the square gradient and + // apply the gradient transpose to v. + Eigen::SimplicialLLT> solver; + Eigen::SparseMatrix squareGradient = gradient.transpose() * gradient; + solver.compute(squareGradient); + Eigen::VectorXd g = solver.solve(gradient.transpose() * v); + + // Find the iso value and shift g by it. + Eigen::MatrixXd sigmaProduct = fullWeightMatrix * g; + printf("Sigma product size: %d, %d\n", sigmaProduct.rows(), sigmaProduct.cols()); + double sigma = sigmaProduct.sum() / n; + printf("Sigma: %f\n", sigma); + g.array() -= 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 + //////////////////////////////////////////////////////////////////////////// + igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F); + + printf("Done\n"); }