Jiayi(Eris) Zhang HW1

src/fd_grad.cpp

@@ -1,4 +1,9 @@
 #include "fd_grad.h"
+#include "igl/cat.h"
+#include "fd_partial_derivative.h"
+#include <iostream>
+
+using namespace Eigen;
 
 void fd_grad(
   const int nx,
@@ -7,7 +12,10 @@ void fd_grad(
   const double h,
   Eigen::SparseMatrix<double> & G)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  SparseMatrix<double> Dx, Dy, Dz, buffer;
+  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);
+  igl::cat(1, Dx, Dy, buffer);
+  igl::cat(1, buffer, Dz, G);
 }

src/fd_interpolate.cpp

@@ -1,4 +1,8 @@
 #include "fd_interpolate.h"
+#include <iostream>
+
+using namespace std;
+using namespace Eigen;
 
 void fd_interpolate(
   const int nx,
@@ -9,7 +13,43 @@ void fd_interpolate(
   const Eigen::MatrixXd & P,
   Eigen::SparseMatrix<double> & W)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  vector<Triplet<double>> triplets;
+  triplets.reserve(P.rows() * 8);
+
+  for (int i = 0; i < P.rows(); i++) {
+    RowVector3d diff = P.row(i) - corner;
+
+    double xr = diff(0) * 1.0 / h;
+    double yr = diff(1) * 1.0 / h;
+    double zr = diff(2) * 1.0 / h;
+
+    int xi = floor(xr);
+    double xd = xr - xi;
+
+    int yi = floor(yr);
+    double yd = yr - yi;
+
+    int zi = floor(zr);
+    double zd = zr - zi;
+
+    // xi, yi, zi
+    triplets.push_back(Triplet<double>(i, nx * ny * zi + nx * yi + xi, (1 - xd) * (1 - yd) * (1 - zd)));
+    // xi + 1, yi, zi
+    triplets.push_back(Triplet<double>(i, nx * ny * zi + nx * yi + (xi + 1), xd * (1 - yd) * (1 - zd)));
+    // xi, yi + 1, zi
+    triplets.push_back(Triplet<double>(i, nx * ny * zi + nx * (yi + 1) + xi, (1 - xd) * yd * (1 - zd)));
+    // xi, yi, zi + 1
+    triplets.push_back(Triplet<double>(i, nx * ny * (zi + 1) + nx * yi + xi, (1 - xd) * (1 - yd) * zd));
+    // xi + 1, yi + 1, zi
+    triplets.push_back(Triplet<double>(i, nx * ny * zi + nx * (yi + 1) + (xi + 1), xd * yd * (1 - zd)));
+    // xi + 1, yi, zi + 1
+    triplets.push_back(Triplet<double>(i, nx * ny * (zi + 1) + nx * yi + (xi + 1), xd * (1 - yd) * zd));
+    // xi, yi + 1, zi + 1
+    triplets.push_back(Triplet<double>(i, nx * ny * (zi + 1) + nx * (yi + 1) + xi, (1 - xd) * yd * zd));
+    // xi + 1, yi + 1, zi + 1
+    triplets.push_back(Triplet<double>(i, nx * ny * (zi + 1) + nx * (yi + 1) + (xi + 1), xd * yd * zd));
+  }
+
+  W.resize(P.rows(), nx*ny*nz);
+  W.setFromTriplets(triplets.begin(), triplets.end());
 }

src/fd_partial_derivative.cpp

@@ -1,5 +1,8 @@
 #include "fd_partial_derivative.h"
 
+using namespace Eigen;
+using namespace std;
+
 void fd_partial_derivative(
   const int nx,
   const int ny,
@@ -8,7 +11,43 @@ void fd_partial_derivative(
   const int dir,
   Eigen::SparseMatrix<double> & D)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  int nxd = nx;
+  int nyd = ny;
+  int nzd = nz;
+
+  if (dir == 0) {
+    nxd--;
+  }
+  else if (dir == 1) {
+    nyd--;
+  }
+  else if (dir == 2) {
+    nzd--;
+  }
+
+  vector<Triplet<double>> triplets;
+  triplets.reserve(D.rows()*2);
+
+  for (int i = 0; i < nxd; i++) {
+    for (int j = 0; j < nyd; j++) {
+      for (int k = 0; k < nzd; k++) {
+        // same for all the three cases
+        int cur = nxd * nyd * k + nxd * j + i;
+        triplets.push_back(Triplet<double>(cur, nx * ny * k + nx * j + i, - 1.0));
+
+        if (dir == 0) {
+          triplets.push_back(Triplet<double>(cur, nx * ny * k + nx * j + (i + 1), 1.0));
+        }
+        else if (dir == 1) {
+          triplets.push_back(Triplet<double>(cur, nx * ny * k + nx * (j + 1) + i, 1.0));
+        }
+        else if (dir == 2) {
+          triplets.push_back(Triplet<double>(cur, nx * ny * (k + 1) + nx * j + i, 1.0));
+        }
+      }
+    }
+  }
+
+  D.resize(nxd*nyd*nzd, nx*ny*nz);
+  D.setFromTriplets(triplets.begin(), triplets.end());
 }

src/poisson_surface_reconstruction.cpp

@@ -1,12 +1,18 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include "fd_grad.h"
+#include "fd_interpolate.h"
+#include <iostream>
+
+using namespace std;
+using namespace Eigen;
 
 void poisson_surface_reconstruction(
-    const Eigen::MatrixXd & P,
-    const Eigen::MatrixXd & N,
-    Eigen::MatrixXd & V,
-    Eigen::MatrixXi & F)
+    const MatrixXd & P,
+    const MatrixXd & N,
+    MatrixXd & V,
+    MatrixXi & F)
 {
   ////////////////////////////////////////////////////////////////////////////
   // Construct FD grid, CONGRATULATIONS! You get this for free!
@@ -23,30 +29,54 @@ void poisson_surface_reconstruction(
   // 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.);
+  RowVector3d corner = P.colwise().minCoeff().array()-pad*h;
+  // Grid dimensions should be at least 3
+  nx = max((P.col(0).maxCoeff()-P.col(0).minCoeff()+(2.*pad)*h)/h,3.);
+  ny = max((P.col(1).maxCoeff()-P.col(1).minCoeff()+(2.*pad)*h)/h,3.);
+  nz = 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++) 
+  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);
+         x.row(ind) = corner + h * RowVector3d(i,j,k);
       }
     }
   }
-  Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
+  VectorXd g = VectorXd::Zero(nx*ny*nz);
+
+  VectorXd vx, vy, vz, v;
+  SparseMatrix<double> Wx, Wy, Wz;
+  //notice the shift of corner here
+  fd_interpolate(nx - 1, ny, nz, h, corner + RowVector3d(0.5 * h, 0, 0), P, Wx);
+  vx = Wx.transpose() * N.col(0);
+  fd_interpolate(nx, ny - 1, nz, h, corner + RowVector3d(0, 0.5 * h, 0), P, Wy);
+  vy = Wy.transpose() * N.col(1);
+  fd_interpolate(nx, ny, nz - 1, h, corner + RowVector3d(0, 0, 0.5 * h), P, Wz);
+  vz = Wz.transpose() * N.col(2);
+
+  v.resize(vx.rows() + vy.rows() + vz.rows()); // remember to initialize size of v here
+  v << vx, vy, vz;
+
+  SparseMatrix<double> G, left, W; // product of sparsematrices are also sparse
+  VectorXd right;
+  fd_grad(nx, ny, nz, h, G);
+  left = G.transpose() * G;
+  right = G.transpose() * v;
+  BiCGSTAB<SparseMatrix<double>> solver;
+  solver.compute(left);
+  g = solver.solve(right);
+
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+  // compute sigma and then shift g
+  g = (g.array() - (MatrixXd::Ones(1, n) * W * g).value() * 1.0 / n).matrix();
 
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
 
   ////////////////////////////////////////////////////////////////////////////
   // Run black box algorithm to compute mesh from implicit function: this