Kenan Deng, HW01

src/fd_grad.cpp

@@ -1,4 +1,7 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <iostream>
+#include <igl/cat.h>
 
 void fd_grad(
   const int nx,
@@ -9,5 +12,18 @@ void fd_grad(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+    //std::cout<<"Grad Here!"<<std::endl;
+    Eigen::SparseMatrix<double> Dx( (nx-1) * ny * nz, nx * ny * nz);
+    Eigen::SparseMatrix<double> Dy( nx * (ny-1) * nz, nx * ny * nz);
+    Eigen::SparseMatrix<double> Dz( nx * ny * (nz-1), nx * ny * nz);
+    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);
+
+    Eigen::SparseMatrix<double> M(Dx.rows() + Dy.rows(), Dx.cols());
+    igl::cat(1, Dx, Dy, M);
+    igl::cat(1, M, Dz, G);
+    //std::cout<<"Grad Done!"<<std::endl;
   ////////////////////////////////////////////////////////////////////////////
 }
+

src/fd_interpolate.cpp

@@ -1,5 +1,7 @@
 #include "fd_interpolate.h"
-
+#include <cmath>
+#include <vector>
+#include <iostream>
 void fd_interpolate(
   const int nx,
   const int ny,
@@ -11,5 +13,32 @@ void fd_interpolate(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+  W.resize(P.rows(), nx * ny * nz);
+  std::vector<Eigen::Triplet<double>> tripletList;
+  tripletList.reserve(8 * P.rows());
+  for (int i = 0; i< P.rows(); i++){
+    double x = std::floor((P(i, 0) - corner(0,0))/h); 
+    double y = std::floor((P(i, 1) - corner(0,1))/h);
+    double z = std::floor((P(i, 2) - corner(0,2))/h);
+    double dx = (P(i, 0) - corner(0,0))/h - x;
+    double dy = (P(i, 1) - corner(0,1))/h - y;
+    double dz = (P(i, 2) - corner(0,2))/h - z;
+
+
+    tripletList.push_back(Eigen::Triplet<double> (i, x + y *nx + z * nx*ny, (1-dx) * (1-dy) * (1-dz)));
+    tripletList.push_back(Eigen::Triplet<double> (i, (x+1) + y *nx + z * nx*ny, dx * (1-dy) * (1-dz)));
+
+    tripletList.push_back(Eigen::Triplet<double> (i, x + (y+1) *nx + z * nx*ny, (1-dx) * dy * (1-dz)));
+    tripletList.push_back(Eigen::Triplet<double> (i, x + y * nx + (z+1) * nx*ny, (1-dx) * (1-dy) * dz));
+
+    tripletList.push_back(Eigen::Triplet<double> (i, (x+1) + (y+1) *nx + z * nx*ny, dx * dy * (1-dz)));
+    tripletList.push_back(Eigen::Triplet<double> (i, (x+1) + y *nx + (z+1) * nx*ny, dx  * (1-dy) * dz));
+
+    tripletList.push_back(Eigen::Triplet<double> (i, x + (y+1) *nx + (z+1) * nx*ny, (1-dx) * dy * dz));
+    tripletList.push_back(Eigen::Triplet<double> (i, (x+1) + (y+1) * nx + (z+1) * nx*ny, dx * dy * dz));
+
+  }
   ////////////////////////////////////////////////////////////////////////////
+
+  W.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/fd_partial_derivative.cpp

@@ -1,5 +1,5 @@
 #include "fd_partial_derivative.h"
-
+#include <iostream>
 void fd_partial_derivative(
   const int nx,
   const int ny,
@@ -10,5 +10,45 @@ void fd_partial_derivative(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+    //std::cout<<"Here!"<<std::endl;
+    if (dir == 0) {
+        D.resize((nx-1)*ny*nz, nx*ny*nz);
+        for (int i = 0; i < (nx-1); ++i){
+            for (int j = 0; j <ny; ++j){
+                for (int k = 0; k<nz; ++k){
+                    D.insert(i + j * (nx-1) + k * (nx-1)  * ny, i + j * nx + k * nx * ny) = -1;
+                    D.insert(i + j * (nx-1) + k * (nx-1)  * ny, i+1 + j * nx + k * nx * ny) = 1;
+                }
+
+            }
+
+        }
+    } else if (dir == 1){
+        D.resize(nx*(ny-1)*nz, 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){
+                    D.insert(i + j * nx + k * nx  * (ny-1), i + j * nx + k * nx * ny) = -1;
+                    D.insert(i + j * nx + k * nx  * (ny-1), i + (j+1) * nx + k * nx * ny) = 1;
+                }
+
+            }
+
+        }
+    } else if (dir == 2){
+        D.resize(nx*ny*(nz-1), 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){
+                    D.insert(i + j * nx + k * nx  * ny, i + j * nx + k * nx * ny) = -1;
+                    D.insert(i + j * nx + k * nx  * ny, i + j * nx + (k+1) * nx * ny) = 1;
+                }
+
+            }
+
+        }
+    }
+    //std::cout<<"Done!"<<std::endl;
+
   ////////////////////////////////////////////////////////////////////////////
 }

src/poisson_surface_reconstruction.cpp

@@ -1,7 +1,9 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
-
+#include "fd_grad.h"
+#include "fd_interpolate.h"
+#include <iostream>
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
     const Eigen::MatrixXd & N,
@@ -43,14 +45,49 @@ void poisson_surface_reconstruction(
     }
   }
   Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
-
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+
+  //Distribute normal 
+  Eigen::VectorXd v = Eigen::VectorXd::Zero((nx-1) * ny * nz + nx * (ny-1) * nz + nx * ny * (nz -1));
+  //Vx
+  Eigen::SparseMatrix<double> W_x(n, (nx-1) * ny * nz);
+  Eigen::VectorXd v_x((nx-1) * ny * nz);
+  fd_interpolate(nx-1,ny,nz,h,corner + h * 0.5 * Eigen::RowVector3d(1,0,0), P,W_x);
+  v_x = W_x.transpose() * N.col(0);
+  //Vy
+  Eigen::SparseMatrix<double> W_y(n, nx * (ny-1) * nz);
+  Eigen::VectorXd v_y(nx * (ny-1) * nz);
+  fd_interpolate(nx,ny-1,nz,h,corner + h * 0.5 * Eigen::RowVector3d(0,1,0), P,W_y);
+  v_y = W_y.transpose() * N.col(1);
+  //Vz
+  Eigen::SparseMatrix<double> W_z(n, nx * ny * (nz-1));
+  Eigen::VectorXd v_z(nx * ny * (nz-1));
+  fd_interpolate(nx,ny,nz-1,h,corner + h * 0.5 * Eigen::RowVector3d(0,0,1), P,W_z);
+  v_z = W_z.transpose() * N.col(2);
+  v << v_x, v_y, v_z;
+  //solve
+  Eigen::SparseMatrix<double> G((nx-1) * ny * nz + nx * (ny-1) * nz + nx * ny * (nz -1), 
+    nx * ny * nz);
+  fd_grad(nx, ny, nz, h, G);
+
+  //Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> cg;
+  Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> cg;
+  cg.compute(G.transpose() * G);
+  Eigen::VectorXd b(nx * ny * nz);
+  b = G.transpose() * v;
+  g = cg.solve(b);
+  //sigma
+  Eigen::SparseMatrix<double> W_sigma(n, nx * ny * nz);
+  fd_interpolate(nx, ny, nz, h, corner, P, W_sigma);
+  double sigma = 1.00/n * (Eigen::MatrixXd::Ones(1, n) * W_sigma *g).value();
+  g = g - sigma * Eigen::MatrixXd::Ones(nx*ny*nz,1);
   ////////////////////////////////////////////////////////////////////////////
 
   ////////////////////////////////////////////////////////////////////////////
   // Run black box algorithm to compute mesh from implicit function: this
   // function always extracts g=0, so "pre-shift" your g values by -sigma
   ////////////////////////////////////////////////////////////////////////////
+
   igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F);
 }