Minhua Zhu HW1

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,13 @@ void fd_grad(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+  Eigen::SparseMatrix<double> Dx, Dy, Dz, Dxy;
+
+  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,Dxy);
+  igl::cat(1,Dxy,Dz,G);
+
 }

src/fd_interpolate.cpp

@@ -1,4 +1,6 @@
 #include "fd_interpolate.h"
+#include <cmath>
+#include <iostream>
 
 void fd_interpolate(
   const int nx,
@@ -12,4 +14,40 @@ void fd_interpolate(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+
+  W.resize(P.rows(),nx * ny * nz);
+  typedef Eigen::Triplet<double> T;
+  std::vector< T > tripletList;
+  tripletList.reserve(8*P.rows());
+
+  for (int i = 0; i < P.rows(); i++) {
+    // distance
+    double dx = P(i,0) - corner[0];
+    double dy = P(i,1) - corner[1];
+    double dz = P(i,2) - corner[2];
+
+    // position of grids
+    int gx = floor(dx / h);
+    int gy = floor(dy / h);
+    int gz = floor(dz / h);
+
+    // relative distance
+    double rx = (dx - gx*h) / h;
+    double ry = (dy - gy*h) / h ;
+    double rz = (dz - gz*h) / h;
+
+    // eight point of one cube
+    tripletList.push_back(T(i, gx+nx*(gy+gz*ny), (1-rx)*(1-ry)*(1-rz)));
+    tripletList.push_back(T(i, gx+1+nx*(gy+gz*ny), rx*(1-ry)*(1-rz)));
+    tripletList.push_back(T(i, gx+nx*(gy+1+gz*ny), (1-rx)*ry*(1-rz)));
+
+    tripletList.push_back(T(i, gx+nx*(gy+(gz+1)*ny), (1-rx)*(1-ry)*rz));
+    tripletList.push_back(T(i, gx+1+nx*(gy+1+gz*ny), rx*ry*(1-rz)));
+    tripletList.push_back(T(i, gx+1+nx*(gy+(gz+1)*ny), rx*(1-ry)*rz));
+
+    tripletList.push_back(T(i, gx+nx*(gy+1+(gz+1)*ny), (1-rx)*ry*rz));
+    tripletList.push_back(T(i, gx+1+nx*(gy+1+(gz+1)*ny), rx*ry*rz));
+
+  }
+  W.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/fd_partial_derivative.cpp

@@ -11,4 +11,55 @@ void fd_partial_derivative(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
-}
+  int m = 0;
+  if (dir == 0) { // nx-1
+    int change_x = nx - 1;
+    m = change_x*ny*nz;
+    D.resize(m, nx*ny*nz);
+    for (int i=0; i<change_x;i++){
+      for (int j=0; j<ny; j++){
+        for (int k=0; k<nz; k++){
+          int index = i + nx * (j + ny * k);
+          int dindex = i + change_x * (j + ny * k); // change nx
+          int nindex = i + 1 + nx * (j + ny * k); // change i
+          D.coeffRef(dindex,index) = - 1/h;
+          D.coeffRef(dindex,nindex) = 1/h;
+        }
+      }
+    }
+  }
+
+  else if (dir == 1)
+  { // ny - 1
+    m = nx*(ny-1)*nz;
+    D.resize(m, nx*ny*nz);
+    for (int i=0; i<nx;i++){
+      for (int j=0; j<ny-1; j++){
+        for (int k=0; k<nz; k++){
+          int index = i + nx * (j + ny * k);
+          int dindex = i + nx * (j + (ny-1) * k); // change ny
+          int nindex = i + nx * (j + 1 + ny * k); // change j
+          D.coeffRef(dindex,index) = - 1/h;
+          D.coeffRef(dindex,nindex) = 1/h;
+        }
+      }
+    }
+  }
+
+  else { // dir == 2 nz-1
+    m = nx*ny*(nz-1);
+    D.resize(m, nx*ny*nz);
+    for (int i=0; i<nx;i++){
+      for (int j=0; j<ny; j++){
+        for (int k=0; k<nz-1; k++){
+          int index = i + nx * (j + ny * k);
+          int dindex = i + nx * (j + ny * k); // change nz
+          int nindex = i + nx * (j + ny * (k+1)); // change z
+          D.coeffRef(dindex,index) = - 1/h;
+          D.coeffRef(dindex,nindex) = 1/h;
+        }
+      }
+    }
+  }
+
+ }

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,12 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include "fd_partial_derivative.h"
+#include "fd_grad.h"
+#include "fd_interpolate.h" 
+#include <Eigen/Sparse>
+#include <Eigen/IterativeLinearSolvers>
+#include <iostream>
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -44,6 +50,8 @@ void poisson_surface_reconstruction(
   }
   Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
 
+
+
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
@@ -52,5 +60,45 @@ void poisson_surface_reconstruction(
   // Run black box algorithm to compute mesh from implicit function: this
   // function always extracts g=0, so "pre-shift" your g values by -sigma
   ////////////////////////////////////////////////////////////////////////////
+
+  Eigen::SparseMatrix<double> Weightx(n,(nx-1)*ny*nz);
+  Eigen::RowVector3d cornerx = corner;
+  cornerx[0] = corner[0] + h/2;
+  fd_interpolate(nx-1,ny,nz,h,cornerx,P,Weightx);
+
+  Eigen::SparseMatrix<double> Weighty(n,nx*(ny-1)*nz);
+  Eigen::RowVector3d cornery = corner;
+  cornery[1] = corner[1] + h/2;
+  fd_interpolate(nx,ny-1,nz,h,corner,P,Weighty);
+
+  Eigen::SparseMatrix<double> Weightz(n,nx*ny*(nz-1));
+  Eigen::RowVector3d cornerz = corner;
+  cornerz[2] = corner[2] + h/2;
+  fd_interpolate(nx,ny,nz-1,h,corner,P,Weightz);
+
+  Eigen::SparseMatrix<double> W(n,nx*ny*nz);
+  fd_interpolate(nx,ny,nz,h,corner,P, W);
+
+  Eigen::SparseMatrix<double> Gradient;
+  fd_grad(nx,ny,nz,h,Gradient);
+
+  Eigen::VectorXd vx(n,1), vy(n,1), vz(n,1);
+  vx=Weightx.transpose()*N.col(0);
+  vy=Weighty.transpose()*N.col(1);
+  vz=Weightz.transpose()*N.col(2);
+  Eigen::VectorXd B((nx-1)*ny*nz+ nx*(ny-1)*nz+ nx*ny*(nz-1));
+  B << vx,vy,vz;
+
+  Eigen::SparseMatrix<double> A = Gradient.transpose()*Gradient;
+  Eigen::VectorXd b = Gradient.transpose()*B;
+  //Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> bicg;
+  Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::UpLoType::Lower | Eigen::UpLoType::Upper> bicg;
+  g = bicg.compute(A).solve(b);
+
+  Eigen::VectorXd i(g.size());
+  i.setOnes();
+  double sigma = (1/(n*1.0)) * (W*g).sum();
+  g = g - sigma*i;
+
   igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F);
 }