Ziheng Liang HW01

include/fd_grad.h

@@ -19,4 +19,10 @@ void fd_grad(
   const int nz,
   const double h,
   Eigen::SparseMatrix<double> & G);
+
+void stack_sparse_matric(
+  Eigen::SparseMatrix<double> &A,
+  Eigen::SparseMatrix<double> &B,
+  Eigen::SparseMatrix<double> &C, 
+  Eigen::SparseMatrix<double> &out);
 #endif

src/fd_grad.cpp

@@ -1,4 +1,6 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <iostream>
 
 void fd_grad(
   const int nx,
@@ -9,5 +11,43 @@ void fd_grad(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+  int dx = (nx-1)*ny*nz;
+  int dy = nx*(ny-1)*nz;
+  int dz = nx*ny*(nz-1);
+
+  Eigen::SparseMatrix<double> Dx = Eigen::SparseMatrix<double>(dx, nx*ny*nz);
+  Eigen::SparseMatrix<double> Dy = Eigen::SparseMatrix<double>(dy, nx*ny*nz);
+  Eigen::SparseMatrix<double> Dz = Eigen::SparseMatrix<double>(dz, nx*ny*nz);
+
+  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);
+  stack_sparse_matric(Dx, Dy, Dz, G);
   ////////////////////////////////////////////////////////////////////////////
 }
+
+//concatenate method copy from stackoverflow
+void stack_sparse_matric(
+  Eigen::SparseMatrix<double> &A,
+  Eigen::SparseMatrix<double> &B,
+  Eigen::SparseMatrix<double> &C, 
+  Eigen::SparseMatrix<double> &out) 
+{
+  std::vector<Eigen::Triplet<double> > tripletList;
+  for (int k = 0; k < A.outerSize(); ++k) {
+    for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it) {
+      tripletList.push_back(Eigen::Triplet<double>(it.row(), it.col(), it.value()));
+    }
+  }
+  for (int k = 0; k < B.outerSize(); ++k) {
+    for (Eigen::SparseMatrix<double>::InnerIterator it(B, k); it; ++it) {
+      tripletList.push_back(Eigen::Triplet<double>(A.rows() + it.row(), it.col(), it.value()));
+    }
+  }
+  for (int k = 0; k < C.outerSize(); ++k) {
+    for (Eigen::SparseMatrix<double>::InnerIterator it(C, k); it; ++it) {
+      tripletList.push_back(Eigen::Triplet<double>(A.rows() + B.rows() + it.row(), it.col(), it.value()));
+    }
+  }
+  out.setFromTriplets(tripletList.begin(), tripletList.end());
+}

src/fd_interpolate.cpp

@@ -1,4 +1,7 @@
 #include "fd_interpolate.h"
+#include <iostream>
+#include <cmath>
+#include <string>
 
 void fd_interpolate(
   const int nx,
@@ -11,5 +14,62 @@ void fd_interpolate(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+  int i_int, j_int, k_int;
+  double i, j, k;
+  double w_x, w_y, w_z;
+  // P -= corner;
+  // Eigen::MatrixXd temp = ((P).rowwise() - corner) / h;
+  // Eigen::MatrixXi P_int = temp.cast <int> ();
+  // Eigen::MatrixXd P_frac = temp - P_int.cast<double> ();
+
+  std::vector<Eigen::Triplet<double> > tripletList;
+
+  for (int row=0; row<P.rows(); row++) {
+
+    i = (P(row,0) - corner(0)) / h;
+    j = (P(row,1) - corner(1)) / h;
+    k = (P(row,2) - corner(2)) / h;
+
+    i_int = std::floor(i);
+    j_int = std::floor(j);
+    k_int = std::floor(k);
+
+    w_x = i - i_int;
+    w_y = j - j_int;
+    w_z = k - k_int;
+
+    // const auto col = i_int + nx*(j_int + k_int * ny);
+    // W.insert(row, col) = (w_x*w_y*w_z);
+    // W.insert(row, col+1) = ((1-w_x)*w_y*w_z);
+    // W.insert(row, col+nx) = (w_x*(1-w_y)*w_z);
+    // W.insert(row, col+nx*ny) = (w_x*w_y*(1-w_z));
+    // W.insert(row, col+1+nx) = ((1-w_x)*(1-w_y)*w_z);
+    // W.insert(row, col+1+nx*ny) = ((1-w_x)*w_y*(1-w_z));
+    // W.insert(row, col+nx+nx*ny) = (w_x*(1-w_y)*(1-w_z));
+    // W.insert(row, col+1+nx+nx*ny) = ((1-w_x)*(1-w_y)*(1-w_z));
+
+    const auto col = i_int + nx*(j_int + k_int * ny);
+
+    // tripletList.push_back(Eigen::Triplet<double>(row, col, P_frac(row,0) * P_frac(row,1) * P_frac(row,2)));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+1, (1-P_frac(row,0)) * P_frac(row,1) * P_frac(row,2)));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+nx, P_frac(row,0) * (1-P_frac(row,1)) * P_frac(row,2)));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+nx*ny, P_frac(row,0) * P_frac(row,1) * (1-P_frac(row,2))));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+1+nx, (1-P_frac(row,0)) * (1-P_frac(row,1)) * P_frac(row,2)));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+1+nx*ny, (1-P_frac(row,0)) * P_frac(row,1) * (1-P_frac(row,2))));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+nx+nx*ny, (P_frac(row,0) * (1-P_frac(row,1)) * (1-P_frac(row,2)))));
+    // tripletList.push_back(Eigen::Triplet<double>(row, col+1+nx+nx*ny, (1-P_frac(row,0)) * (1-P_frac(row,1)) * (1-P_frac(row,2))));
+
+
+    tripletList.push_back(Eigen::Triplet<double>(row, col, (1-w_x) * (1-w_y) * (1-w_z)));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+1, w_x * (1-w_y) * (1-w_z)));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+nx, (1-w_x) * w_y * (1-w_z)));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+nx*ny, (1-w_x) * (1-w_y) * w_z));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+1+nx, w_x * w_y * (1-w_z)));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+1+nx*ny, w_x * (1-w_y) * w_z));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+nx+nx*ny, (1-w_x) * w_y * w_z));
+    tripletList.push_back(Eigen::Triplet<double>(row, col+1+nx+nx*ny, w_x * w_y * w_z));
+
+  }
+  W.setFromTriplets(tripletList.begin(), tripletList.end());
   ////////////////////////////////////////////////////////////////////////////
 }

src/fd_partial_derivative.cpp

@@ -10,5 +10,32 @@ void fd_partial_derivative(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+  int x_diff = 0, y_diff = 0, z_diff = 0;
+  if (dir == 0) {
+  	x_diff = 1;
+  }
+  else if (dir == 1) {
+  	y_diff = 1;
+  }
+  else if (dir == 2) {
+  	z_diff = 1;
+  }
+
+  D = Eigen::SparseMatrix<double>((nx-x_diff)*(ny-y_diff)*(nz-z_diff), nx*ny*nz);
+  const auto offset = (x_diff + nx*(y_diff + z_diff * ny));
+  for(int i = 0; i < nx - x_diff; i++) 
+  {
+    for(int j = 0; j < ny - y_diff; j++)
+    {
+      for(int k = 0; k < nz - z_diff; k++)
+      {
+        // Convert subscript to index
+        const auto stagger_idx = i + j * (nx-x_diff) + k * (ny-y_diff) * (nx-x_diff);
+        const auto original_idx = i + j * nx + k * ny * nx;
+        D.insert(stagger_idx, original_idx) = -1/h;
+        D.insert(stagger_idx, original_idx + offset) = 1/h;
+      }
+    }
+  }
   ////////////////////////////////////////////////////////////////////////////
 }

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,10 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include <Eigen/Sparse>
+#include "fd_grad.h"
+#include "fd_interpolate.h"
+#include <iostream>
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -47,6 +51,37 @@ void poisson_surface_reconstruction(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+  int dx = (nx-1)*ny*nz;
+  int dy = nx*(ny-1)*nz;
+  int dz = nx*ny*(nz-1);
+
+  Eigen::SparseMatrix<double> Wx(P.rows(), dx);
+  Eigen::SparseMatrix<double> Wy(P.rows(), dy);
+  Eigen::SparseMatrix<double> Wz(P.rows(), dz);
+  Eigen::SparseMatrix<double> W(P.rows(), nx*ny*nz);
+  Eigen::SparseMatrix<double> G(dx+dy+dz, nx*ny*nz);
+
+  fd_grad(nx, ny, nz, h, G);
+  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);
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+  Eigen::VectorXd vx, vy, vz;
+  vx = Wx.transpose() * N.col(0);
+  vy = Wy.transpose() * N.col(1);
+  vz = Wz.transpose() * N.col(2);
+
+  Eigen::VectorXd v(dx+dy+dz);
+  v << vx, vy, vz;
+  Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
+  Eigen::SparseMatrix<double> G_square = G.transpose() * G;
+  Eigen::VectorXd Gv = G.transpose() * v;
+  solver.compute(G_square);
+  g = solver.solve(Gv);
+
+  double sigma = (W * g).sum() / n;
+  g = (g.array() - sigma).matrix();
 
   ////////////////////////////////////////////////////////////////////////////
   // Run black box algorithm to compute mesh from implicit function: this