From 1232acda78f973e692ea60ab994bf7cfd098b265 Mon Sep 17 00:00:00 2001 From: David Acuna Date: Wed, 26 Sep 2018 22:32:14 -0400 Subject: [PATCH] hw 01.. --- src/fd_grad.cpp | 28 +++-- src/fd_interpolate.cpp | 72 +++++++++++-- src/fd_partial_derivative.cpp | 61 +++++++++-- src/poisson_surface_reconstruction.cpp | 135 ++++++++++++++++--------- 4 files changed, 217 insertions(+), 79 deletions(-) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..1c52bee 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,13 +1,23 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" +#include "igl/cat.h" 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, Dy, Dz, buff; + 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); + + //matrix: G = [Dx;Dy;Dz] + + igl::cat(1, Dx, Dy, buff); + igl::cat(1, buff, Dz, G); + } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..8d8e962 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,15 +1,65 @@ #include "fd_interpolate.h" +#include +#include "igl/floor.h" + +//aux function that helps to keep the code clean and easy to follow. +int fast_indexing(Eigen::RowVector3d &p_floor, int i, int j, int k, int nx, int ny) { + return (p_floor(0) + i) + (nx * (p_floor(1) + j)) + (nx * ny * (p_floor(2) + k)); + +} 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) { + + + typedef Eigen::Triplet et; + std::vector TList; + TList.reserve(8 * P.rows()); + + + for (int i = 0; i < P.rows(); i++) { + // Let's normalize all the points. + Eigen::RowVector3d p = (1. / h) * (P.row(i) - corner); + Eigen::RowVector3d p_floor; + + // Let's compute p_0 + igl::floor(p, p_floor); + + Eigen::RowVector3d diff = p - p_floor; + + //This could be done in a for-loop but it obscures the code. + // In any case, it is just 8 lines. + TList.push_back(et(i, fast_indexing(p_floor, 0, 0, 0, nx, ny), + (1 - diff[0]) * (1 - diff[1]) * (1 - diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 0, 0, 1, nx, ny), + (1 - diff[0]) * (1 - diff[1]) * (diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 0, 1, 0, nx, ny), + (1 - diff[0]) * (diff[1]) * (1 - diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 0, 1, 1, nx, ny), + (1 - diff[0]) * (diff[1]) * (diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 1, 0, 0, nx, ny), + (diff[0]) * (1 - diff[1]) * (1 - diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 1, 0, 1, nx, ny), + (diff[0]) * (1 - diff[1]) * (diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 1, 1, 0, nx, ny), + (diff[0]) * (diff[1]) * (1 - diff[2]))); + + TList.push_back(et(i, fast_indexing(p_floor, 1, 1, 1, nx, ny), + (diff[0]) * (diff[1]) * (diff[2]))); + } + + W.setFromTriplets(TList.begin(), TList.end()); + } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..2ba4792 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,14 +1,55 @@ #include "fd_partial_derivative.h" +#include + +int fast_indexing(int i, int j, int k, int nx, int ny) { + return i + (nx * j) + (nx * ny * k); + +} 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) { + + int nxm = nx - (dir == 0); + int nym = ny - (dir == 1); + int nzm = nz - (dir == 2); + + + // this triplet is + // by using this we can save time while setting the sparseM D + typedef Eigen::Triplet et; + std::vector TList; + + //number of element to insert + TList.reserve(nxm * nym * nzm); + + D.resize(nxm * nym * nzm, nx * ny * nz); + + for (int i = 0; i < nxm; i++) { + for (int j = 0; j < nym; j++) { + for (int k = 0; k < nzm; k++) { + TList.push_back(et(fast_indexing(i, j, k, nxm, nym), fast_indexing(i, j, k, nx, ny), -1)); + et n; + switch (dir) { + case 0: + n = et(fast_indexing(i, j, k, nxm, nym), fast_indexing(i + 1, j, k, nx, ny), 1); + break; + case 1: + n = et(fast_indexing(i, j, k, nxm, nym), fast_indexing(i, j + 1, k, nx, ny), 1); + break; + case 2: + n = et(fast_indexing(i, j, k, nxm, nym), fast_indexing(i, j, k + 1, nx, ny), 1); + break; + + } + TList.push_back(n); + } + } + } + D.setFromTriplets(TList.begin(), TList.end()); + } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..3bf1c30 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,56 +1,93 @@ #include "poisson_surface_reconstruction.h" #include #include +#include "fd_interpolate.h" +#include "fd_grad.h" +#include "igl/cat.h" void poisson_surface_reconstruction( - const Eigen::MatrixXd & P, - const Eigen::MatrixXd & N, - 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); - } + const Eigen::MatrixXd &P, + const Eigen::MatrixXd &N, + 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); + Eigen::VectorXd g = Eigen::VectorXd::Zero(nx * ny * nz); + + Eigen::SparseMatrix Wx(P.rows(), (nx - 1) * ny * nz); + Eigen::SparseMatrix Wy(P.rows(), nx * (ny - 1) * nz); + Eigen::SparseMatrix Wz(P.rows(), nx * ny * (nz - 1)); + + // We need to shift the corners. + fd_interpolate(nx - 1, ny, nz, h, corner + Eigen::RowVector3d(h / 2, 0., 0.), P, Wx); + fd_interpolate(nx, ny - 1, nz, h, corner + Eigen::RowVector3d(0., h / 2., 0.), P, Wy); + fd_interpolate(nx, ny, nz - 1, h, corner + Eigen::RowVector3d(0, 0., h / 2), P, Wz); + + // Computing V + Eigen::VectorXd vx = Wx.transpose() * N.col(0); + Eigen::VectorXd vy = Wy.transpose() * N.col(1); + Eigen::VectorXd vz = Wz.transpose() * N.col(2); + + //Merging v + Eigen::VectorXd buff, v; + igl::cat(1, vx, vy, buff); + igl::cat(1, buff, vz, v); + + // Computing G + Eigen::SparseMatrix G; + fd_grad(nx, ny, nz, h, G); + + // Solving the system + std::cout << "Solving the system...this may take a bit.." << std::endl; + Eigen::BiCGSTAB> solver; + solver.compute(G.transpose() * G); + g = solver.solve(G.transpose() * v); + std::cout << "[+] Solver Error:" << solver.error() << std::endl; + + //Computing sigma + Eigen::SparseMatrix W(P.rows(), nx * ny * nz); + fd_interpolate(nx, ny, nz, h, corner, P, W); + + //according to the paper: sigma = 1/n*1^T*W*g => \R^(1xn)*\R^(n,nx*ny*nz)*\R^(nx*ny*nz,1) + //This is equivalent to take the mean between \R^(n,nx*ny*nz)*\R^(nx*ny*nz,1)=mean(\R^(n,1)) + double sigma = (W * g).mean(); + + //pre-shifting g values by subtracting sigma + g = g - Eigen::VectorXd::Ones(g.size()) * sigma; + + + //////////////////////////////////s////////////////////////////////////////// + // 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); }