Hw1 Jienan Yao Submission

src/fd_grad.cpp

@@ -1,13 +1,38 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
 
 void fd_grad(
-  const int nx,
-  const int ny,
-  const int nz,
-  const double h,
-  Eigen::SparseMatrix<double> & G)
+	const int nx,
+	const int ny,
+	const int nz,
+	const double h,
+	Eigen::SparseMatrix<double> & G)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+	////////////////////////////////////////////////////////////////////////////
+	// Add your code here
+	Eigen::SparseMatrix<double> Gx, Gy, Gz;
+
+	// Calculate partial derivative matrix for each x, y, z direction
+	fd_partial_derivative(nx, ny, nz, h, 0, Gx);
+	fd_partial_derivative(nx, ny, nz, h, 1, Gy);
+	fd_partial_derivative(nx, ny, nz, h, 2, Gz);
+
+	// Concatenate the matrix vertically
+	Eigen::SparseMatrix<double> M(Gx.rows() + Gy.rows() + Gz.rows(), Gx.cols());
+	M.reserve(Gx.nonZeros() + Gy.nonZeros() + Gz.nonZeros());
+	for (int c = 0; c < Gx.cols(); c++) {
+		M.startVec(c);
+		for (Eigen::SparseMatrix<double>::InnerIterator itx(Gx, c); itx; ++itx) {
+			M.insertBack(itx.row(), c) = itx.value();
+		}
+		for (Eigen::SparseMatrix<double>::InnerIterator ity(Gy, c); ity; ++ity) {
+			M.insertBack(ity.row() + Gx.rows(), c) = ity.value();
+		}
+		for (Eigen::SparseMatrix<double>::InnerIterator itz(Gz, c); itz; ++itz) {
+			M.insertBack(itz.row() + Gx.rows() + Gy.rows(), c) = itz.value();
+		}		
+	}
+	////////////////////////////////////////////////////////////////////////////
+	G = M;
+	M.finalize();
 }

src/fd_interpolate.cpp

@@ -1,4 +1,5 @@
 #include "fd_interpolate.h"
+#include <math.h>
 
 void fd_interpolate(
   const int nx,
@@ -11,5 +12,44 @@ void fd_interpolate(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+	W.setZero();
+	std::vector<Eigen::Triplet<double>> triplets;
+	triplets.reserve(8 * P.rows());
+
+	for (int r = 0; r < P.rows(); r++) {
+		double x = P(r, 0);
+		double y = P(r, 1);
+		double z = P(r, 2);
+
+		int i, j, k;
+		// x-component
+		i = floor((x - corner(0)) / h);
+		j = floor((y - corner(1)) / h);
+		k = floor((z - corner(2)) / h);
+
+		// The coordinates of the left-bottom-front
+		// point of the grid surrounding P
+		double grid_corner_x, grid_corner_y, grid_corner_z;
+		grid_corner_x = corner(0) + h * i;
+		grid_corner_y = corner(1) + h * j;
+		grid_corner_z = corner(2) + h * k;
+
+		// The fraction of the coordinate of P relative to the left-bottom-front
+		// point of the surrounding grid
+		double tx, ty, tz;
+		tx = (x - grid_corner_x) / h;
+		ty = (y - grid_corner_y) / h;
+		tz = (z - grid_corner_z) / h;
+		triplets.push_back({r, i + j * nx + k * nx * ny, (1 - tx) * (1 - ty) * (1 - tz)});
+		triplets.push_back({r, i + (j + 1) * nx + k * ny * nx, (1 - tx) * ty * (1 - tz)});
+		triplets.push_back({r, i + j * nx + (k + 1) * ny * nx, (1 - tx) * (1 - ty) * (tz)});
+		triplets.push_back({r, i + (j + 1) * nx + (k + 1) * ny * nx, (1 - tx) * (ty) * (tz)});
+		triplets.push_back({r, i + 1 + (j)* nx + (k)* ny * nx, tx * (1 - ty) * (1 - tz)});
+		triplets.push_back({r, i + 1 + (j + 1) * nx + (k)* ny * nx, (tx) * (ty) * (1 - tz)});
+		triplets.push_back({r, i + 1 + (j)* nx + (k + 1) * ny * nx, (tx) * (1 - ty) * (tz)});
+		triplets.push_back({r, i + 1 + (j + 1) * nx + (k + 1) * ny * nx, (tx) * (ty) * (tz)});
+	}
+	W.setFromTriplets(triplets.begin(), triplets.end());
+
   ////////////////////////////////////////////////////////////////////////////
 }

src/fd_partial_derivative.cpp

@@ -10,5 +10,75 @@ void fd_partial_derivative(
 {
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
+
+	int m;
+	if (dir == 0) {
+		m = (nx - 1)*ny*nz;
+	}
+	else if (dir == 1) {
+		m = nx * (ny - 1)*nz;
+	}
+	else {
+		m = nx * ny*(nz - 1);
+	}
+	Eigen::SparseMatrix<double> Dx(m, nx*ny*nz);
+	Dx.setZero();
+	std::vector<Eigen::Triplet<double>> triplets;
+	triplets.reserve(m*2);
+
+	// Form the finite difference matrices based
+	int curr_row = 0;
+	int curr_col = 0;
+	if (dir == 0) {
+		for (int k = 0; k < nz; k++) {
+			for (int j = 0; j < ny; j++) {
+				for (int i = 0; i < nx; i++) {
+					if (i == nx - 1) {
+						curr_col += 1;
+					}
+					else {
+						triplets.push_back({ curr_row, curr_col  , -1  });
+						triplets.push_back({ curr_row, curr_col + 1,  1  });
+						curr_row++; curr_col++;
+					}
+				}
+			}
+		}
+	}
+	else if (dir == 1) {
+		for (int k = 0; k < nz; k++) {
+			for (int j = 0; j < ny; j++) {
+				for (int i = 0; i < nx; i++) {
+					if (j == ny - 1) {
+						curr_col += 1;
+					}	else {
+						triplets.push_back({ curr_row, curr_col  , -1 });
+						triplets.push_back({ curr_row, curr_col + nx,  1  });
+						curr_row++; curr_col++;
+					}
+				}
+			}
+		}
+	}
+	else {
+		for (int k = 0; k < nz; k++) {
+			for (int j = 0; j < ny; j++) {
+				for (int i = 0; i < nx; i++) {
+					if (k == nz - 1) {
+						curr_col += 1;
+					}
+					else {
+
+						triplets.push_back({ curr_row, i + j * nx + k * nx*ny  , -1  });
+						triplets.push_back({ curr_row, i + j * nx + k * nx*ny + nx * ny,  1  });
+						curr_row++;
+					}
+				}
+			}
+		}
+	}
+
+	Dx.setFromTriplets(triplets.begin(), triplets.end());
+	D = Dx;
   ////////////////////////////////////////////////////////////////////////////
 }

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,8 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include "fd_grad.h"
+#include "fd_interpolate.h"
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -47,6 +49,52 @@ void poisson_surface_reconstruction(
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////
+  //Interpolation matrices on the staggered grids and the primary grid
+  Eigen::SparseMatrix<double> Wx(P.rows(), (nx - 1) * ny * nz);
+  Eigen::SparseMatrix<double> Wy(P.rows(), (ny - 1) * nx * nz);
+  Eigen::SparseMatrix<double> Wz(P.rows(), (nz - 1) * nx * ny);
+  Eigen::SparseMatrix<double> W(P.rows(), nx*ny*nz);
+  Eigen::SparseMatrix<double> G;
+  // Corner points of the staggered grids
+  Eigen::RowVector3d cornerx, cornery, cornerz;
+  cornerx(0) = corner(0) + h/2;  cornerx(1) = corner(1);  cornerx(2) = corner(2);
+  cornery(0) = corner(0);  cornery(1) = corner(1) + h/2;  cornery(2) = corner(2);
+  cornerz(0) = corner(0);  cornerz(1) = corner(1);  cornerz(2) = corner(2) + h/2;
+
+  // Calculate interpolation matrix on the staggered grids
+  fd_interpolate(nx - 1, ny, nz, h, cornerx, P, Wx);
+  fd_interpolate(nx, ny - 1, nz, h, cornery, P, Wy);
+  fd_interpolate(nx, ny, nz - 1, h, cornerz, P, Wz);
+
+  // Calculate interpolation matrix W on the primary grid
+  fd_interpolate(nx, ny, nz, h, corner, P, W);
+  Eigen::MatrixXd vx(P.rows(), 1), vy(P.rows(), 1), vz(P.rows(), 1);
+  vx = Wx.transpose() * N.col(0);
+  vy = Wy.transpose() * N.col(1);
+  vz = Wz.transpose() * N.col(2);
+
+  // Concatenate vx, vy and vz
+  Eigen::VectorXd v(vx.rows()+vy.rows()+vz.rows());
+  for (int i = 0; i < vx.rows(); i++) {
+	  v(i) = vx(i);
+  }
+  for (int i = 0; i < vy.rows(); i++) {
+	  v(i+vx.rows()) = vy(i);
+  }
+  for (int i = 0; i < vz.rows(); i++) {
+	  v(i+vx.rows()+vy.rows()) = vz(i);
+  }
+
+  fd_grad(nx, ny, nz, h, G);
+
+  Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower | Eigen::Upper> solver;
+  g = solver.compute(G.transpose() * G).solve(G.transpose() * v);
+
+  double sigma = (1/(double) P.rows()) * (W * g).sum();
+
+  Eigen::VectorXd ones(g.size());
+  ones.setOnes();
+  g = g - sigma * ones;
 
   ////////////////////////////////////////////////////////////////////////////
   // Run black box algorithm to compute mesh from implicit function: this