[Meng Wei] hw01

src/fd_grad.cpp

@@ -1,4 +1,7 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <igl/cat.h>
+#include <iostream>
 
 void fd_grad(
   const int nx,
@@ -7,7 +10,16 @@ void fd_grad(
   const double h,
   Eigen::SparseMatrix<double> & G)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  // Construct gradient:
+  Eigen::SparseMatrix<double> Dx, Dy, Dz, Inter;
+  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);
+
+  // Cat:
+  Inter.resize(Dx.rows() + Dy.rows(), Dx.cols());
+  igl::cat(1, Dx, Dy, Inter);
+  G.resize(Inter.rows() + Dz.rows(), Dx.cols());
+  igl::cat(1, Inter, Dz, G);
+  G.makeCompressed();
 }

src/fd_interpolate.cpp

@@ -1,4 +1,8 @@
 #include "fd_interpolate.h"
+#include <iostream>
+
+// Idea Inspired by Eigen library manual of setFromTriplets() method
+typedef Eigen::Triplet<double> T;
 
 void fd_interpolate(
   const int nx,
@@ -9,7 +13,65 @@ void fd_interpolate(
   const Eigen::MatrixXd & P,
   Eigen::SparseMatrix<double> & W)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  // Construct a matrix(P.rows, P.cols) with all rows equal to corner
+  Eigen::MatrixXd corner_matrix(P.rows(), P.cols());
+  corner_matrix << Eigen::MatrixXd::Ones(P.rows(), P.cols());
+
+  corner_matrix.col(0) << corner_matrix.col(0) * corner.col(0);
+  corner_matrix.col(1) << corner_matrix.col(1) * corner.col(1);
+  corner_matrix.col(2) << corner_matrix.col(2) * corner.col(2);
+
+  // Compute relative position in the grid
+  Eigen::MatrixXd position(P.rows(), P.cols());
+  position << (P - corner_matrix) / h;
+
+  Eigen::MatrixXd grid(P.rows(), P.cols());
+  grid << (position.cast <int> ()).cast <double> ();
+
+  Eigen::MatrixXd abso(P.rows(), P.cols());
+  abso << position - grid;
+
+
+  // Idea inspired by Eigen library manual of setFromTriplets() method
+  std::vector<T> triplet_list; 
+  triplet_list.reserve(P.rows() * 8);
+
+  for (int i = 0; i < P.rows(); i ++) {
+    // x, y, z
+    triplet_list.push_back(
+      T(i, grid(i, 0) + grid(i, 1)*nx + grid(i, 2)*nx*ny, (1-abso(i, 0))*(1-abso(i, 1))*(1-abso(i, 2)) )
+    );
+    // x+1, y, z
+    triplet_list.push_back(
+      T(i, grid(i, 0)+1 + grid(i, 1)*nx + grid(i, 2)*nx*ny, abso(i, 0)*(1-abso(i, 1))*(1-abso(i, 2)))
+    );
+    // x, y+1, z
+    triplet_list.push_back(
+      T(i, grid(i, 0) + (grid(i, 1)+1)*nx + grid(i, 2)*nx*ny,  (1-abso(i, 0))*abso(i,1)*(1-abso(i, 2)))
+    );
+    // x+1, y+1, z
+    triplet_list.push_back(
+      T(i, grid(i, 0)+1 + (grid(i, 1)+1)*nx + grid(i, 2)*nx*ny, abso(i, 0)*abso(i, 1)*(1-abso(i, 2)))
+    );
+    // x, y, z+1
+    triplet_list.push_back(
+      T(i, grid(i, 0) + grid(i, 1)*nx + (grid(i, 2)+1)*nx*ny, (1-abso(i, 0))*(1-abso(i, 1))*abso(i, 2) )
+    );
+    // x+1, y, z+1
+    triplet_list.push_back(
+      T(i, grid(i, 0)+1 + grid(i, 1)*nx + (grid(i, 2)+1)*nx*ny, abso(i, 0)*(1-abso(i, 1))*abso(i, 2))
+    );
+    // x, y+1, z+1
+    triplet_list.push_back(
+      T(i, grid(i, 0) + (grid(i, 1)+1)*nx + (grid(i, 2)+1)*nx*ny, (1-abso(i, 0))*abso(i, 1)*abso(i, 2))
+    );
+    // x+1, y+1, z+1
+    triplet_list.push_back(
+      T(i, grid(i, 0)+1 + (grid(i, 1)+1)*nx + (grid(i, 2)+1)*nx*ny, abso(i, 0)*abso(i, 1)*abso(i, 2))
+    );
+
+  }
+  W.resize(P.rows(), nx*ny*nz);
+  W.setFromTriplets(triplet_list.begin(), triplet_list.end());
+  W.makeCompressed();
 }

src/fd_partial_derivative.cpp

@@ -1,4 +1,8 @@
 #include "fd_partial_derivative.h"
+#include <iostream>
+
+// Idea Inspired by Eigen library manual of setFromTriplets() method
+typedef Eigen::Triplet<double> T;
 
 void fd_partial_derivative(
   const int nx,
@@ -8,7 +12,77 @@ void fd_partial_derivative(
   const int dir,
   Eigen::SparseMatrix<double> & D)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  // Prepare new shape:
+  int new_x = (dir == 0) ? nx - 1 : nx;
+  int new_y = (dir == 1) ? ny - 1 : ny;
+  int new_z = (dir == 2) ? nz - 1 : nz;
+
+  int num_point = nx * ny * nz;
+  int total = new_x * new_y * new_z;
+
+  // Idea inspired by Eigen library manual of setFromTriplets() method
+  std::vector<T> triplet_list; 
+
+  // Each row contains two values (-1, 1):
+  triplet_list.reserve(total*2);
+  switch (dir) {
+    case 0:
+      for(int k = 0; k < new_z; k ++) {
+        for(int j = 0; j < new_y; j ++) {
+          for(int i = 0; i < new_x; i ++) {
+            triplet_list.push_back(
+              T(i + new_x*j + new_x*new_y*k,
+                // skip previous nodes:
+                i + nx*j + nx*ny*k, 
+                -1/h));
+            triplet_list.push_back(
+              T(i + new_x*j + new_x*new_y*k,
+                // move 1 column:
+                i+1 + nx*j + nx*ny*k,
+                1/h));
+          }
+        }
+      }
+      break;
+    case 1:
+      for(int k = 0; k < new_z; k ++) {
+        for(int j = 0; j < new_y; j ++) {
+          for(int i = 0; i < new_x; i ++) {
+            triplet_list.push_back(
+              T(i + new_x*j + new_x*new_y*k,
+                // skip previous nodes:
+                i + nx*j + nx*ny*k, 
+                -1/h));
+            triplet_list.push_back(
+              T(i + new_x*j + new_x*new_y*k,
+                // move 1 row:
+                i + nx*(j+1) + nx*ny*k,
+                1/h));
+          }
+        }
+      }
+      break;
+    case 2:
+      for(int k = 0; k < new_z; k ++) {
+        for(int j = 0; j < new_y; j ++) {
+          for(int i = 0; i < new_x; i ++) {
+            triplet_list.push_back(
+              T(i + new_x*j + new_x*new_y*k,
+                // skip previous nodes:
+                i + nx*j + nx*ny*k, 
+                -1/h));
+            triplet_list.push_back(
+              T(i + new_x*j + new_x*new_y*k,
+                // move 1 plane:
+                i + nx*j + nx*ny*(k+1),
+                1/h));
+          }
+        }
+      }
+      break;
+  }
+
+  D.resize(total, num_point);
+  D.setFromTriplets(triplet_list.begin(), triplet_list.end());
+  D.makeCompressed();
 }

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,10 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include <iostream>
+#include "fd_interpolate.h"
+#include "fd_partial_derivative.h"
+#include "fd_grad.h"
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -45,12 +49,46 @@ void poisson_surface_reconstruction(
   Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
 
   ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
+  // distribute the given normals N 
+  Eigen::SparseMatrix<double> Wx, Wy, Wz;
+
+  // Moving P to the left by h/2 is the same as moving corner to the right by h/2
+  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);
+
+  Eigen::VectorXd Vx(Wx.cols());
+  Vx << Wx.transpose() * N.col(0);
+  Eigen::VectorXd Vy(Wy.cols());
+  Vy << Wy.transpose() * N.col(1);
+  Eigen::VectorXd Vz(Wz.cols());
+  Vz << Wz.transpose() * N.col(2);
+  Eigen::VectorXd my_vec(Vx.size() + Vy.size() + Vz.size());
+  my_vec << Vx, Vy, Vz;
+
+
+  // Compute Gradient:
+  Eigen::SparseMatrix<double> G;
+  fd_grad(nx, ny, nz, h, G);
+
+
+  // Compute g: G.T G g = G.T v
+  // Idea inspired by BiCGSTAB library
+  Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper> cg;
+  cg.compute(G.transpose() * G);
+  g = cg.solve(G.transpose() * my_vec);
+
   ////////////////////////////////////////////////////////////////////////////
 
   ////////////////////////////////////////////////////////////////////////////
-  // Run black box algorithm to compute mesh from implicit function: this
-  // function always extracts g=0, so "pre-shift" your g values by -sigma
+  // Find W on the primary (non-staggered) grid:
+  Eigen::SparseMatrix<double> W;
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+  // Compute sigma:
+  double sigma = (W * g).sum() / W.rows();
+  g = g - Eigen::VectorXd::Ones(g.size()) * sigma;
+
   ////////////////////////////////////////////////////////////////////////////
   igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F);
 }