HW01 David Acuna

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<double> & G)
-{
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+        const int nx,
+        const int ny,
+        const int nz,
+        const double h,
+        Eigen::SparseMatrix<double> &G) {
+
+
+    Eigen::SparseMatrix<double> 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);
+
 }

src/fd_interpolate.cpp

@@ -1,15 +1,65 @@
 #include "fd_interpolate.h"
+#include <math.h>
+#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<double> & 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<double> &W) {
+
+
+    typedef Eigen::Triplet<double> et;
+    std::vector<et> 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());
+
 }

src/fd_partial_derivative.cpp

@@ -1,14 +1,55 @@
 #include "fd_partial_derivative.h"
+#include <iostream>
+
+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<double> & D)
-{
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+        const int nx,
+        const int ny,
+        const int nz,
+        const double h,
+        const int dir,
+        Eigen::SparseMatrix<double> &D) {
+
+    int nxm = nx - (dir == 0);
+    int nym = ny - (dir == 1);
+    int nzm = nz - (dir == 2);
+
+
+    // this triplet is <row,col,value>
+    // by using this we can save time while setting the sparseM D
+    typedef Eigen::Triplet<double> et;
+    std::vector<et> 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());
+
 }

src/poisson_surface_reconstruction.cpp

@@ -1,56 +1,93 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#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<double> Wx(P.rows(), (nx - 1) * ny * nz);
+    Eigen::SparseMatrix<double> Wy(P.rows(), nx * (ny - 1) * nz);
+    Eigen::SparseMatrix<double> 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<double> 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<Eigen::SparseMatrix<double>> solver;
+    solver.compute(G.transpose() * G);
+    g = solver.solve(G.transpose() * v);
+    std::cout << "[+] Solver Error:" << solver.error() << std::endl;
+
+    //Computing sigma
+    Eigen::SparseMatrix<double> 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);
 }