Lawson Fulton

include/poisson_surface_reconstruction.h

@@ -1,7 +1,7 @@
 #ifndef POISSON_SURFACE_RECONSTRUCTION_H
 #define POISSON_SURFACE_RECONSTRUCTION_H
 #include <Eigen/Core>
-
+#include <Eigen/Sparse>
 // Takes input sample points P and input normals N 
 // and gives a watertight mesh using a simplified
 // version of [Kazhdan et. al 2006]

src/fd_grad.cpp

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

src/fd_interpolate.cpp

@@ -1,5 +1,10 @@
 #include "fd_interpolate.h"
 
+#include <iostream>
+#include <vector>
+
+#include "igl/floor.h"
+
 void fd_interpolate(
   const int nx,
   const int ny,
@@ -9,7 +14,42 @@ void fd_interpolate(
   const Eigen::MatrixXd & P,
   Eigen::SparseMatrix<double> & W)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  // Get the grid corner indices for each input point
+  Eigen::MatrixXd P_rel_grid = (P.rowwise() - corner) / h;
+
+  Eigen::MatrixXi P_corners_ijk;
+  igl::floor(P_rel_grid, P_corners_ijk);
+
+  // These are the first three weights
+  Eigen::MatrixXd P_rel_corner = P_rel_grid - P_corners_ijk.cast<double>();
+
+  auto ind = [&nx, &ny](int i, int j, int k) {
+    return i + nx * (j + k * ny);
+  };
+
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> tripletList;
+  tripletList.reserve(8 * P.rows());
+  for(int l = 0; l < P.rows(); l++) 
+  {
+    int i = P_corners_ijk(l, 0);
+    int j = P_corners_ijk(l, 1);
+    int k = P_corners_ijk(l, 2);
+
+    double wx = P_rel_corner(l, 0);
+    double wy = P_rel_corner(l, 1);
+    double wz = P_rel_corner(l, 2);
+
+    tripletList.push_back(T(l, ind(i+0, j+0, k+0), (1. - wx) * (1. - wy) * (1. - wz)));
+    tripletList.push_back(T(l, ind(i+0, j+0, k+1), (1. - wx) * (1. - wy) * wz));
+    tripletList.push_back(T(l, ind(i+0, j+1, k+0), (1. - wx) * wy        * (1. - wz)));
+    tripletList.push_back(T(l, ind(i+0, j+1, k+1), (1. - wx) * wy        * wz));
+    tripletList.push_back(T(l, ind(i+1, j+0, k+0), wx        * (1. - wy) * (1. - wz)));
+    tripletList.push_back(T(l, ind(i+1, j+0, k+1), wx        * (1. - wy) * wz));
+    tripletList.push_back(T(l, ind(i+1, j+1, k+0), wx        * wy        * (1. - wz)));
+    tripletList.push_back(T(l, ind(i+1, j+1, k+1), wx        * wy        * wz));
+  }
+
+  W.resize(P.rows(), nx*ny*nz);
+  W.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/fd_partial_derivative.cpp

@@ -8,7 +8,43 @@ void fd_partial_derivative(
   const int dir,
   Eigen::SparseMatrix<double> & D)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  int is_x = dir == 0;
+  int is_y = dir == 1;
+  int is_z = dir == 2;
+
+  int nxd = nx - 1 * is_x;
+  int nyd = ny - 1 * is_y;
+  int nzd = nz - 1 * is_z;
+
+  auto ind = [&nx, &ny](int i, int j, int k) {
+    return i + nx * (j + k * ny);
+  };
+
+  auto ind_d = [&nxd, &nyd](int i, int j, int k) {
+    return i + nxd * (j + k * nyd);
+  };
+
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> tripletList;
+  tripletList.reserve(nxd*nyd*nzd);
+
+  for(int i = 0; i < nxd; i++) 
+  {
+    for(int j = 0; j < nyd; j++)
+    {
+      for(int k = 0; k < nzd; k++)
+      {
+        int ip = i + 1 * is_x;
+        int jp = j + 1 * is_y;
+        int kp = k + 1 * is_z;
+
+        int row = ind_d(i, j, k);
+        tripletList.push_back(T(row, ind(i,j,k),     1./h));
+        tripletList.push_back(T(row, ind(ip,jp,kp), -1./h));
+      }
+    }
+  }
+
+  D.resize(nxd*nyd*nzd, nx*ny*nz);
+  D.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,9 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include <iostream>
+#include "fd_interpolate.h"
+#include "fd_grad.h"
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -47,7 +50,41 @@ void poisson_surface_reconstruction(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+  Eigen::RowVector3d dx(h/2., 0., 0.);
+  Eigen::RowVector3d dy(0., h/2., 0.);
+  Eigen::RowVector3d dz(0., 0., h/2.);
+  Eigen::SparseMatrix<double> Wx;
+  fd_interpolate(nx - 1, ny, nz, h, corner + dx, P.rowwise() + dx, Wx);
+  Eigen::VectorXd vx = Wx.transpose() * N.col(0);
+  std::cout << vx.size() << std::endl;
 
+  Eigen::SparseMatrix<double> Wy;
+  fd_interpolate(nx, ny - 1, nz, h, corner + dy, P.rowwise() + dy, Wy);
+  Eigen::VectorXd vy = Wy.transpose() * N.col(1);
+
+  Eigen::SparseMatrix<double> Wz;
+  fd_interpolate(nx, ny, nz - 1, h, corner + dz, P.rowwise() + dz, Wz);
+  Eigen::VectorXd vz = Wz.transpose() * N.col(2);
+
+  Eigen::VectorXd v(vx.size() + vy.size() + vz.size());
+  v << vx, vy, vz;
+
+  Eigen::SparseMatrix<double> G;
+  fd_grad(nx, ny, nz, h, G);
+  std::cout << "G rows: " << G.rows() << " G cols: " << G.cols() << std::endl;
+
+  std::cout << "Solving..." << std::endl;
+  Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
+  solver.compute(G.transpose() * G);
+  g = solver.solve(G.transpose() * v);
+  std::cout << "#iterations:     " << solver.iterations() << std::endl;
+  std::cout << "estimated error: " << solver.error()      << std::endl;
+
+  Eigen::SparseMatrix<double> W;
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+  double sigma = (W * g).mean();
+  g = -(g - Eigen::VectorXd::Ones(g.size()) * 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