HW1 - John (Junrui) Xu

src/fd_grad.cpp

@@ -1,4 +1,6 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <igl/cat.h>
 
 void fd_grad(
   const int nx,
@@ -10,4 +12,16 @@ void fd_grad(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+
+    Eigen::SparseMatrix<double> Dx;
+    Eigen::SparseMatrix<double> Dy;
+    Eigen::SparseMatrix<double> Dz;
+
+    fd_partial_derivative(nx, ny, nz, h, 1, Dx);
+    fd_partial_derivative(nx, ny, nz, h, 2, Dy);
+    fd_partial_derivative(nx, ny, nz, h, 3, Dz);
+
+    Eigen::SparseMatrix<double> t;
+    igl::cat(1, Dx, Dy, t);
+    igl::cat(1, t, Dz, G);
 }

src/fd_interpolate.cpp

@@ -12,4 +12,33 @@ void fd_interpolate(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+
+    W.resize(P.rows(), nx * ny * nz);
+
+    std::vector<Eigen::Triplet<double>> buffer;
+
+    for(int i = 0; i < P.rows(); i++){
+        int x0 = (P(i, 0) - corner(0)) / h;
+        int x1 = x0 + 1;
+        double dx = (P(i, 0) - (x0 * h + corner(0))) / h;
+
+        int y0 = (P(i, 1) - corner(1)) / h;
+        int y1 = y0 + 1;
+        double dy = (P(i, 1) - (y0 * h + corner(1))) / h;
+
+        int z0 = (P(i, 2) - corner(2)) / h;
+        int z1 = z0 + 1;
+        double dz = (P(i, 2) - (z0 * h + corner(2))) / h;
+
+        buffer.push_back(Eigen::Triplet<double>(i, x0 + y0*nx + z0*ny*nx,(1.0 - dx)*(1.0 - dy)*(1.0 - dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x1 + y0*nx + z0*ny*nx,(      dx)*(1.0 - dy)*(1.0 - dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x0 + y1*nx + z0*ny*nx,(1.0 - dx)*(      dy)*(1.0 - dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x1 + y1*nx + z0*ny*nx,(      dx)*(      dy)*(1.0 - dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x0 + y0*nx + z1*ny*nx,(1.0 - dx)*(1.0 - dy)*(      dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x1 + y0*nx + z1*ny*nx,(      dx)*(1.0 - dy)*(      dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x0 + y1*nx + z1*ny*nx,(1.0 - dx)*(      dy)*(      dz)));
+        buffer.push_back(Eigen::Triplet<double>(i, x1 + y1*nx + z1*ny*nx,(      dx)*(      dy)*(      dz)));
+    }
+
+    W.setFromTriplets(buffer.begin(), buffer.end());
 }

src/fd_partial_derivative.cpp

@@ -1,4 +1,5 @@
 #include "fd_partial_derivative.h"
+#include <assert.h>
 
 void fd_partial_derivative(
   const int nx,
@@ -11,4 +12,48 @@ void fd_partial_derivative(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+
+    int nx_D = nx, ny_D = ny, nz_D = nz;
+
+    switch (dir) {
+    case 1:
+        nx_D--;
+        break;
+    case 2:
+        ny_D--;
+        break;
+    case 3:
+        nz_D--;
+        break;
+    default:
+        assert(0);
+    }
+
+    D.resize(nx_D * ny_D * nz_D, nx * ny * nz);
+
+    for (int i = 0; i < nx_D; i++) {
+        for (int j = 0; j < ny_D; j++) {
+            for (int k = 0; k < nz_D; k++) {
+
+                int staggered_ijk = i + j * nx_D + k * ny_D * nx_D;
+                int l_prev = i + j * nx + k * ny * nx;
+                int l_curr;
+
+                switch (dir) {
+                case 1:
+                    l_curr = (i + 1) + j * nx + k * ny * nx;
+                    break;
+                case 2:
+                    l_curr = i + (j + 1) * nx + k * ny * nx;
+                    break;
+                case 3:
+                    l_curr = i + j * nx + (k + 1) * ny * nx;
+                    break;
+                }
+
+                D.insert(staggered_ijk, l_prev) = -1;
+                D.insert(staggered_ijk, l_curr) = 1;
+            }
+        }
+    }
 }

src/poisson_surface_reconstruction.cpp

@@ -1,4 +1,6 @@
 #include "poisson_surface_reconstruction.h"
+#include "fd_grad.h"
+#include "fd_interpolate.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
 
@@ -48,6 +50,42 @@ void poisson_surface_reconstruction(
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
 
+  // Staggered interpolation matrices
+
+  Eigen::SparseMatrix<double> Wx;
+  Eigen::SparseMatrix<double> Wy;
+  Eigen::SparseMatrix<double> Wz;
+
+  fd_interpolate(nx - 1, ny, nz, h, corner, P, Wx);
+  fd_interpolate(nx, ny - 1, nz, h, corner, P, Wy);
+  fd_interpolate(nx, ny, nz - 1, h, corner, P, Wz);
+
+  // Construct 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);
+  Eigen::VectorXd v(vx.rows() + vy.rows() + vz.rows());
+  v << vx, vy, vz;
+
+  // Construct and solve the linear system
+
+  Eigen::SparseMatrix<double> G;
+  fd_grad(nx, ny, nz, h, G);
+
+  Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver(G.transpose() * G);
+  g = solver.solve(G.transpose() * v);
+
+  // Determine the iso-level sigma
+
+  Eigen::SparseMatrix<double> W;
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+  double sigma = Eigen::VectorXd::Constant(P.rows(), 1.0).transpose()*W*g;
+  sigma /= P.rows();
+
+  g -= Eigen::VectorXd::Constant(g.rows(), sigma);
+
   ////////////////////////////////////////////////////////////////////////////
   // Run black box algorithm to compute mesh from implicit function: this
   // function always extracts g=0, so "pre-shift" your g values by -sigma