Zhiyuan Cheng's HW01

src/fd_grad.cpp

@@ -1,13 +1,19 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <igl/cat.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
-  ////////////////////////////////////////////////////////////////////////////
+	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,15 +1,44 @@
 #include "fd_interpolate.h"
+#include <math.h>
 
 void fd_interpolate(
-  const int nx,
-  const int ny,
-  const int nz,
-  const double h,
-  const Eigen::RowVector3d & corner,
-  const Eigen::MatrixXd & P,
-  Eigen::SparseMatrix<double> & W)
+	const int nx,
+	const int ny,
+	const int nz,
+	const double h,
+	const Eigen::RowVector3d & corner,
+	const Eigen::MatrixXd & P,
+	Eigen::SparseMatrix<double> & W)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
-}
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	tripletList.reserve(P.rows() * 8);
+	for(int i = 0; i < P.rows(); i++) {
+		// relative position to corner
+		double x = (P(i, 0) - corner[0])/h;
+		double y = (P(i, 1) - corner[1])/h;
+		double z = (P(i, 2) - corner[2])/h;
+
+		// previous on grid position
+		int x1 = floor(x);
+		int y1 = floor(y);
+		int z1 = floor(z);
+
+		// next on grid position
+		int x2 = x1 + 1;
+		int y2 = y1 + 1;
+		int z2 = z1 + 1;
+
+		// calculate weight for each of eight surrounding point
+		tripletList.push_back(T(i, (x1 + y1 * nx + z1 * ny * nx), (x2 - x) * (y2 - y) * (z2 - z)));
+		tripletList.push_back(T(i, (x1 + y1 * nx + z2 * ny * nx), (x2 - x) * (y2 - y) * (z - z1)));
+		tripletList.push_back(T(i, (x1 + y2 * nx + z1 * ny * nx), (x2 - x) * (y - y1) * (z2 - z)));
+		tripletList.push_back(T(i, (x1 + y2 * nx + z2 * ny * nx), (x2 - x) * (y - y1) * (z - z1)));
+		tripletList.push_back(T(i, (x2 + y1 * nx + z1 * ny * nx), (x - x1) * (y2 - y) * (z2 - z)));
+		tripletList.push_back(T(i, (x2 + y1 * nx + z2 * ny * nx), (x - x1) * (y2 - y) * (z - z1)));
+		tripletList.push_back(T(i, (x2 + y2 * nx + z1 * ny * nx), (x - x1) * (y - y1) * (z2 - z)));
+		tripletList.push_back(T(i, (x2 + y2 * nx + z2 * ny * nx), (x - x1) * (y - y1) * (z - z1)));
+	}
+	W.resize(P.rows(), nx * ny * nz);
+	W.setFromTriplets(tripletList.begin(), tripletList.end());
+}

src/fd_partial_derivative.cpp

@@ -1,14 +1,54 @@
 #include "fd_partial_derivative.h"
 
 void fd_partial_derivative(
-  const int nx,
-  const int ny,
-  const int nz,
-  const double h,
-  const int dir,
-  Eigen::SparseMatrix<double> & D)
+	const int nx,
+	const int ny,
+	const int nz,
+	const double h,
+	const int dir,
+	Eigen::SparseMatrix<double> & D)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+	// compute number of entries
+	int x = nx;
+	int y = ny;
+	int z = nz;
+	switch (dir){
+		case 0:
+			x--;
+			break;
+		case 1:
+			y--;
+			break;
+		case 2:
+			z--;
+			break;
+	}
+
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	tripletList.reserve(x * y * z * 2);
+
+	for(int i = 0; i < x; i++) {
+		for (int j = 0; j < y; j++) {
+			for (int k = 0; k < z; k++) {
+				int s = i + j * x + k * y * x;
+				int t = i + j * nx + k * ny * nx;
+				tripletList.push_back(T(s, t, -1/h));
+				switch (dir){
+					case 0:
+						t = (i + 1) + j * nx + k * ny * nx;
+						break;
+					case 1:
+						t = i + (j + 1) * nx + k * ny * nx;
+						break;
+					case 2:
+						t = i + j * nx + (k + 1) * ny * nx;
+						break;
+				}
+				tripletList.push_back(T(s, t, 1/h));
+			}
+		}
+	}
+	D.resize(x * y * z, nx * ny * nz);
+	D.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/poisson_surface_reconstruction.cpp

@@ -1,56 +1,91 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include "fd_grad.h"
+#include "fd_interpolate.h"
+#include <igl/cat.h>
 
 void poisson_surface_reconstruction(
-    const Eigen::MatrixXd & P,
-    const Eigen::MatrixXd & N,
-    Eigen::MatrixXd & V,
-    Eigen::MatrixXi & F)
+	const Eigen::MatrixXd & P,
+	const Eigen::MatrixXd & N,
+	Eigen::MatrixXd & V,
+	Eigen::MatrixXi & F)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Construct FD grid, CONGRATULATIONS! You get this for free!
-  ////////////////////////////////////////////////////////////////////////////
-  // number of input points
-  const int n = P.rows();
-  // Grid dimensions
-  int nx, ny, nz;
-  // Maximum extent (side length of bounding box) of points
-  double max_extent =
-    (P.colwise().maxCoeff()-P.colwise().minCoeff()).maxCoeff();
-  // padding: number of cells beyond bounding box of input points
-  const double pad = 8;
-  // choose grid spacing (h) so that shortest side gets 30+2*pad samples
-  double h  = max_extent/double(30+2*pad);
-  // Place bottom-left-front corner of grid at minimum of points minus padding
-  Eigen::RowVector3d corner = P.colwise().minCoeff().array()-pad*h;
-  // Grid dimensions should be at least 3 
-  nx = std::max((P.col(0).maxCoeff()-P.col(0).minCoeff()+(2.*pad)*h)/h,3.);
-  ny = std::max((P.col(1).maxCoeff()-P.col(1).minCoeff()+(2.*pad)*h)/h,3.);
-  nz = std::max((P.col(2).maxCoeff()-P.col(2).minCoeff()+(2.*pad)*h)/h,3.);
-  // Compute positions of grid nodes
-  Eigen::MatrixXd x(nx*ny*nz, 3);
-  for(int i = 0; i < nx; i++) 
-  {
-    for(int j = 0; j < ny; j++)
-    {
-      for(int k = 0; k < nz; k++)
-      {
-         // Convert subscript to index
-         const auto ind = i + nx*(j + k * ny);
-         x.row(ind) = corner + h*Eigen::RowVector3d(i,j,k);
-      }
-    }
-  }
-  Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
+	////////////////////////////////////////////////////////////////////////////
+	// Construct FD grid, CONGRATULATIONS! You get this for free!
+	////////////////////////////////////////////////////////////////////////////
+	// number of input points
+	const int n = P.rows();
+	// Grid dimensions
+	int nx, ny, nz;
+	// Maximum extent (side length of bounding box) of points
+	double max_extent =
+		(P.colwise().maxCoeff() - P.colwise().minCoeff()).maxCoeff();
+	// padding: number of cells beyond bounding box of input points
+	const double pad = 8;
+	// choose grid spacing (h) so that shortest side gets 30+2*pad samples
+	double h = max_extent / double(30 + 2 * pad);
+	// Place bottom-left-front corner of grid at minimum of points minus padding
+	Eigen::RowVector3d corner = P.colwise().minCoeff().array() - pad * h;
+	// Grid dimensions should be at least 3 
+	nx = std::max((P.col(0).maxCoeff() - P.col(0).minCoeff() + (2.*pad)*h) / h, 3.);
+	ny = std::max((P.col(1).maxCoeff() - P.col(1).minCoeff() + (2.*pad)*h) / h, 3.);
+	nz = std::max((P.col(2).maxCoeff() - P.col(2).minCoeff() + (2.*pad)*h) / h, 3.);
+	// Compute positions of grid nodes
+	Eigen::MatrixXd x(nx*ny*nz, 3);
+	for (int i = 0; i < nx; i++)
+	{
+		for (int j = 0; j < ny; j++)
+		{
+			for (int k = 0; k < nz; k++)
+			{
+				// Convert subscript to index
+				const auto ind = i + nx * (j + k * ny);
+				x.row(ind) = corner + h * Eigen::RowVector3d(i, j, k);
+			}
+		}
+	}
+	Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
 
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+	Eigen::SparseMatrix<double> Wx, Wy, Wz, G, W;
+	// get gradient
+	fd_grad(nx, ny, nz, h, G);
 
-  ////////////////////////////////////////////////////////////////////////////
-  // 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);
+	// staggered grid
+	Eigen::RowVector3d corner_x = corner + Eigen::RowVector3d(h / 2, 0, 0);
+	Eigen::RowVector3d corner_y = corner + Eigen::RowVector3d(0, h / 2, 0);
+	Eigen::RowVector3d corner_z = corner + Eigen::RowVector3d(0, 0, h / 2);
+	fd_interpolate((nx - 1), ny, nz, h, corner_x, P, Wx);
+	fd_interpolate(nx, (ny - 1), nz, h, corner_y, P, Wy);
+	fd_interpolate(nx, ny, (nz - 1), h, corner_z, P, Wz);
+
+	// construct v
+	Eigen::VectorXd vx(nx*ny*nz - ny*nz), vy(nx*ny*nz - nx*nz), vz(nx*ny*nz - nx*ny), 
+		vxy(2*nx*ny*nz - ny*nz - nx*nz), v(3*nx*ny*nz - nx*ny - ny*nz - nx*nz);
+	vx = Wx.transpose() * N.col(0);
+	vy = Wy.transpose() * N.col(1);
+	vz = Wz.transpose() * N.col(2);
+	igl::cat(1, vx, vy, vxy);
+	igl::cat(1, vxy, vz, v);
+
+	// solve for g
+	Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver;
+	solver.compute(G.transpose() * G);
+	g = solver.solve(G.transpose() * v);
+
+	// compute iso-value
+	fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+	// shift
+	// g -= (1.0 / n) * Eigen::RowVectorXd::Ones(n) * W * g * Eigen::VectorXd::Ones(nx * ny * nz);
+	double sigma = (W*g).sum()/n;
+  	g = (g.array() - sigma).matrix();
+
+  	////////////////////////////////////////////////////////////////////////////
+	// 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);
+	// std::cout << "V is of size " << V.size() << std::endl;
+	// std::cout << "F is of size " << F.size() << std::endl;
 }