HW01 Submit - Vineeth Bhaskara

src/fd_grad.cpp

@@ -1,4 +1,9 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <iostream>
+
+void matrix_concat_vertical(Eigen::SparseMatrix<double> & G_x, Eigen::SparseMatrix<double> & G_y,
+Eigen::SparseMatrix<double> & G_z, Eigen::SparseMatrix<double> & G);
 
 void fd_grad(
   const int nx,
@@ -9,5 +14,63 @@ void fd_grad(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+
+  Eigen::SparseMatrix<double> G_x;
+  Eigen::SparseMatrix<double> G_y;
+  Eigen::SparseMatrix<double> G_z;
+
+  fd_partial_derivative(nx, ny, nz, h, 0, G_x);
+  fd_partial_derivative(nx, ny, nz, h, 1, G_y);
+  fd_partial_derivative(nx, ny, nz, h, 2, G_z);
+
+  matrix_concat_vertical(G_x, G_y, G_z, G);
+
+
   ////////////////////////////////////////////////////////////////////////////
 }
+
+void matrix_concat_vertical(Eigen::SparseMatrix<double> & G_x, Eigen::SparseMatrix<double> & G_y,
+Eigen::SparseMatrix<double> & G_z, Eigen::SparseMatrix<double> & G) {
+
+// holder for the merged sparse matrices
+  std::vector<Eigen::Triplet<double>> tripletList;
+
+  tripletList.reserve((G_x.rows() + G_y.rows() + G_z.rows())*2);
+
+  // iterate over non-zero elements of G_x
+  for (int k=0; k<G_x.outerSize(); ++k)
+  for (Eigen::SparseMatrix<double>::InnerIterator it(G_x,k); it; ++it)
+  {
+    /*
+      it.value(); // value
+      it.row();   // row index
+      it.col();   // col index (here it is equal to k)
+      it.index(); // inner index, here it is equal to it.row()
+      Ref: https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
+    */
+
+    tripletList.push_back(Eigen::Triplet<double>(it.row(), it.col(), it.value()));
+  }
+
+  // now over G_y
+  for (int k=0; k<G_y.outerSize(); ++k)
+  for (Eigen::SparseMatrix<double>::InnerIterator it(G_y,k); it; ++it)
+  {
+    tripletList.push_back(Eigen::Triplet<double>(G_x.rows() + it.row(), it.col(), it.value()));
+  }
+
+  // now over G_z
+  for (int k=0; k<G_z.outerSize(); ++k)
+  for (Eigen::SparseMatrix<double>::InnerIterator it(G_z,k); it; ++it)
+  {
+
+    tripletList.push_back(Eigen::Triplet<double>(G_x.rows() + G_y.rows() + it.row(), it.col(), it.value()));
+  }
+
+
+  G.resize(G_x.rows()+G_y.rows()+G_z.rows(), G_x.cols());
+  G.setZero();
+
+  G.setFromTriplets(tripletList.begin(), tripletList.end());
+
+}

src/fd_interpolate.cpp

@@ -1,4 +1,9 @@
 #include "fd_interpolate.h"
+#include <cmath>
+#include <vector>
+#include <iostream>
+
+int compute_col_idx_node(int idx_x, int idx_y, int idx_z, int nx, int ny);
 
 void fd_interpolate(
   const int nx,
@@ -11,5 +16,123 @@ void fd_interpolate(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+
+  // Assumption: The Point list P is of dimensions: integer X 3 (That is, coordinates represented by columns).
+
+  // For each of the points, need to find the 1-step grid range in x, y, z to localize the corresponding cube
+  // that is required for trilinear interpolation.
+
+  // requried variables that localize the small unit cube in the grid for a given point 
+  double grid_cube_x_0;
+  double grid_cube_x_1;
+  double grid_cube_y_0;
+  double grid_cube_y_1;
+  double grid_cube_z_0;
+  double grid_cube_z_1;
+
+  // for indexing the grid nodes
+  int grid_cube_x_0_idx;
+  int grid_cube_x_1_idx;
+  int grid_cube_y_0_idx;
+  int grid_cube_y_1_idx;
+  int grid_cube_z_0_idx;
+  int grid_cube_z_1_idx;
+
+  // define intermediate variables
+  double x_d, y_d, z_d;
+  double c_000, c_001, c_100, c_010, c_110, c_101, c_011, c_111;
+
+
+  // the actual matrix
+  std::vector<Eigen::Triplet<double> > tripletList;
+
+  // W would be of dimensions num_points X (nx*ny*nz)
+  // But each row would have utmost 8 non-zero entries
+  // Therefore, reserve the Triplet size to num_points * 8
+  tripletList.reserve(P.rows()*8);
+
+
+
+  for (int i=0; i<P.rows(); i++) {
+    // for each point given, we need the corresponding small cube of the grid to be identified 
+    // as per the question, the corner of the bottom-left-front-most node of the grid is given
+
+    // grid_cube_x_0, grid_cube_x_1 represent the two x-coordinates between which the x-coordinate of the point exists
+    grid_cube_x_0 = corner(0) + double(floor((P(i,0)-corner(0))/h)) * h;
+    grid_cube_x_1 = grid_cube_x_0 + h;
+
+    // similarly for the other coordinates
+    grid_cube_y_0 = corner(1) + double(floor((P(i,1)-corner(1))/h)) * h;
+    grid_cube_y_1 = grid_cube_y_0 + h;
+
+    grid_cube_z_0 = corner(2) + double(floor((P(i,2)-corner(2))/h)) * h;
+    grid_cube_z_1 = grid_cube_z_0 + h;
+
+    // now we have localized the cube to which the point p(i) belong to
+    // let's assign integers to the nodes in the mesh so that it is easy to vectorize as X
+    // this helps in knowing where to fill in the W matrix
+    // Taking Corner index as (0,0,0), we assign indices to all nodes assuming that the indices are all positive (since we need to access it as a matrix index)
+    grid_cube_y_0_idx = int((grid_cube_y_0 - corner(1))/h);
+    grid_cube_y_1_idx = grid_cube_y_0_idx + 1;
+
+    grid_cube_z_0_idx = int((grid_cube_z_0 - corner(2))/h);
+    grid_cube_z_1_idx = grid_cube_z_0_idx + 1;
+
+    grid_cube_x_0_idx = int((grid_cube_x_0 - corner(0))/h);
+    grid_cube_x_1_idx = grid_cube_x_0_idx + 1;
+
+
+
+    // interpolation matrix W's weights are going to based on the distances between the above identified small cube
+    // and the given point - so, W is of dimensions: num_points X (nx * ny * nz)
+    // and given (i, j, k) integer mesh-coordinates of a node, we have it corresponding 
+    // to the column index: i + j * n_x + k * n_y * n_x - in W 
+    // compute the interpolation weights wrt to the cube identified for the point
+
+    // define intermediate variables for computing interpolation
+    x_d = (P(i,0) - grid_cube_x_0)/h;
+    y_d = (P(i,1) - grid_cube_y_0)/h;
+    z_d = (P(i,2) - grid_cube_z_0)/h;
+
+    // compute the weights
+    c_000 = (1.0 - x_d) * (1.0 - y_d) * (1.0 - z_d);
+    c_100 = (x_d) * (1.0 - y_d) * (1.0 - z_d);
+    c_010 = (1.0 - x_d) * (y_d) * (1.0 - z_d);
+    c_110 = (x_d) * (y_d) * (1.0 - z_d);
+    c_001 = (1.0 - x_d) * (1.0 - y_d) * (z_d);
+    c_101 = (x_d) * (1.0 - y_d) * (z_d);
+    c_011 = (1.0 - x_d) * (y_d) * (z_d);
+    c_111 = (x_d) * (y_d) * (z_d);
+
+
+
+    // insert the 8 weights at appropriate columns
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_0_idx, grid_cube_y_0_idx, grid_cube_z_0_idx, nx, ny), c_000));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_1_idx, grid_cube_y_0_idx, grid_cube_z_0_idx, nx, ny), c_100));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_0_idx, grid_cube_y_1_idx, grid_cube_z_0_idx, nx, ny), c_010));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_1_idx, grid_cube_y_1_idx, grid_cube_z_0_idx, nx, ny), c_110));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_0_idx, grid_cube_y_0_idx, grid_cube_z_1_idx, nx, ny), c_001));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_1_idx, grid_cube_y_0_idx, grid_cube_z_1_idx, nx, ny), c_101));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_0_idx, grid_cube_y_1_idx, grid_cube_z_1_idx, nx, ny), c_011));
+    tripletList.push_back(Eigen::Triplet<double>(i, compute_col_idx_node(grid_cube_x_1_idx, grid_cube_y_1_idx, grid_cube_z_1_idx, nx, ny), c_111));
+
+
+    // iterate
+  }
+
+
+
+W.resize(P.rows(), nx*ny*nz);
+W.setZero();
+
+// set the values to the sparse matrices
+W.setFromTriplets(tripletList.begin(), tripletList.end());
+
   ////////////////////////////////////////////////////////////////////////////
 }
+
+// given the mesh node indices, returns the column index in W that corresponds to its interpolation weight contribution
+int compute_col_idx_node(int idx_x, int idx_y, int idx_z, int nx, int ny) {
+  int ind = idx_x + idx_y * nx + idx_z * ny * nx;
+  return ind;
+}

src/fd_partial_derivative.cpp

@@ -1,4 +1,8 @@
 #include "fd_partial_derivative.h"
+#include <vector>
+#include <iostream>
+
+int find_col_idx_node(int idx_x, int idx_y, int idx_z, int nx, int ny);
 
 void fd_partial_derivative(
   const int nx,
@@ -9,6 +13,90 @@ void fd_partial_derivative(
   Eigen::SparseMatrix<double> & D)
 {
   ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
+
+  // the number of rows
+  int num_rows_D;
+
+  // the actual matrix
+  std::vector<Eigen::Triplet<double> > tripletList;
+
+  // some intermediate variables
+  int col_index_before, col_index_after, col_index;
+
+
+  if (dir == 0) {
+    num_rows_D = (nx-1)*ny*nz;
+    // the upper bound for the number is 2 non-zero elements per row
+    tripletList.reserve(num_rows_D*2);
+
+    for (int i=0; i<nx-1; i++) {
+      for (int j=0; j<ny; j++) {
+        for (int k=0; k<nz; k++) {
+          // get the column index for the point and the next point
+          col_index = find_col_idx_node(i, j, k, nx-1, ny);
+          col_index_before = find_col_idx_node(i, j, k, nx, ny);
+          col_index_after = find_col_idx_node(i+1, j, k, nx, ny);
+
+          // assign the derivative
+          tripletList.push_back(Eigen::Triplet<double>(col_index, col_index_before, -1.0));
+          tripletList.push_back(Eigen::Triplet<double>(col_index, col_index_after, +1.0));
+        }
+      }
+    }
+
+
+  } else if (dir == 1) {
+    num_rows_D = nx*(ny-1)*nz;
+    // the upper bound for the number is 2 non-zero elements per row
+    tripletList.reserve(num_rows_D*2);
+
+    for (int i=0; i<nx; i++) {
+      for (int j=0; j<ny-1; j++) {
+        for (int k=0; k<nz; k++) {
+          // get the column index for the point and the next point
+          col_index = find_col_idx_node(i, j, k, nx, ny-1);
+          col_index_before = find_col_idx_node(i, j, k, nx, ny);
+          col_index_after = find_col_idx_node(i, j+1, k, nx, ny);
+
+          // assign the derivative
+          tripletList.push_back(Eigen::Triplet<double>(col_index, col_index_before, -1.0));
+          tripletList.push_back(Eigen::Triplet<double>(col_index, col_index_after, +1.0));
+        }
+      }
+    }
+
+  } else {
+    num_rows_D = nx*ny*(nz-1);
+    // the upper bound for the number is 2 non-zero elements per row
+    tripletList.reserve(num_rows_D*2);
+
+    for (int i=0; i<nx; i++) {
+      for (int j=0; j<ny; j++) {
+        for (int k=0; k<nz-1; k++) {
+          // get the column index for the point and the next point
+          col_index_before = find_col_idx_node(i, j, k, nx, ny);
+          col_index_after = find_col_idx_node(i, j, k+1, nx, ny);
+
+          // assign the derivative
+          tripletList.push_back(Eigen::Triplet<double>(col_index_before, col_index_before, -1.0));
+          tripletList.push_back(Eigen::Triplet<double>(col_index_before, col_index_after, +1.0));
+        }
+      }
+    }
+
+  }
+
+
+
+  D.resize(num_rows_D, nx*ny*nz);
+  D.setZero();
+
+  D.setFromTriplets(tripletList.begin(), tripletList.end());
+
   ////////////////////////////////////////////////////////////////////////////
 }
+
+// given the mesh node indices, returns the corresponding column index
+int find_col_idx_node(int idx_x, int idx_y, int idx_z, int nx, int ny) {
+  return idx_x + idx_y * nx + idx_z * ny * nx;
+}

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,11 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
+#include "fd_interpolate.h"
+#include "fd_grad.h"
 #include <algorithm>
+#include <iostream>
+#include <Eigen/Sparse>
+
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -22,6 +27,12 @@ void poisson_surface_reconstruction(
   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);
+
+  /* Probably, the definition of h should be instead as follows to ensure atleast 30+2*pad samples along the *smallest* dimension:
+        double min_extent = (P.colwise().maxCoeff()-P.colwise().minCoeff()).minCoeff();
+        double h = min_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 
@@ -48,9 +59,61 @@ void poisson_surface_reconstruction(
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
 
+  // to hold the weight matrix for staggered interpolation
+  Eigen::SparseMatrix<double> Wx;
+  Eigen::SparseMatrix<double> Wy;
+  Eigen::SparseMatrix<double> Wz; 
+
+
+  // need to interpolate only the staggered positions
+  fd_interpolate(nx-1, ny, nz, h, corner + (h/2)*Eigen::RowVector3d(1,0,0), P, Wx);
+  fd_interpolate(nx, ny-1, nz, h, corner + (h/2)*Eigen::RowVector3d(0,1,0), P, Wy);
+  fd_interpolate(nx, ny, nz-1, h, corner + (h/2)*Eigen::RowVector3d(0,0,1), P, Wz);
+
+
+  Eigen::MatrixXd N_interpolated_x = Wx.transpose() * N.col(0);
+  Eigen::MatrixXd N_interpolated_y = Wy.transpose() * N.col(1);
+  Eigen::MatrixXd N_interpolated_z = Wz.transpose() * N.col(2);
+
+  // concatenate
+  Eigen::MatrixXd N_interpolated(N_interpolated_x.rows()+N_interpolated_y.rows()+N_interpolated_z.rows(), N_interpolated_x.cols());
+
+  N_interpolated << N_interpolated_x, N_interpolated_y, N_interpolated_z;
+
+
+
+
+  // compute the gradient
+  Eigen::SparseMatrix<double> G;
+
+  fd_grad(nx, ny, nz, h, G);
+
+
+  // solve the system
+  Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
+  solver.compute(G.transpose()*G);
+
+  Eigen::VectorXd g_without_sigma = Eigen::VectorXd::Zero(nx*ny*nz);
+
+  g_without_sigma = solver.solve(G.transpose()*N_interpolated);
+
+  // Shift by the iso-value
+  Eigen::VectorXd n_ones = Eigen::VectorXd::Ones(n);
+
+
+  // to hold the weight matrix for regular grid (g) interpolation
+  Eigen::SparseMatrix<double> W;
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+  double sigma;
+  sigma = (1./n)*n_ones.transpose()*W*g;
+
+  g = g_without_sigma - sigma*Eigen::VectorXd::Ones(nx*ny*nz);
+
   ////////////////////////////////////////////////////////////////////////////
   // 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);
 }
+