Submitting my implementation for this assignment.

include/fd_grad.h

@@ -1,6 +1,7 @@
 #ifndef FD_GRAD_H
 #define FD_GRAD_H
 #include <Eigen/Sparse>
+
 // Construct a gradient matrix for a finite-difference grid
 //
 // Inputs:

main.cpp

@@ -9,6 +9,12 @@
 #include <vector>
 #include <cstdlib>
 
+#ifndef AP_WINDOWS_BUILD
+#define PATH_PREFIX "../"
+#else
+#define PATH_PREFIX ""
+#endif
+
 int main(int argc, char *argv[])
 {
   // Load in points + normals from .pwn file
@@ -18,7 +24,7 @@ int main(int argc, char *argv[])
     std::vector<std::vector<double> > vD;
     std::string line;
     std::fstream in;
-    in.open(argc>1?argv[1]:"../shared/data/hand.pwn");
+    in.open(argc>1?argv[1]:PATH_PREFIX "shared/data/hand.pwn");
     while(in)
     {
       std::getline(in, line);

src/fd_grad.cpp

@@ -1,13 +1,50 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+
+typedef Eigen::Triplet<double> Triplet;
+
+/**
+* @brief Concatenates vertically.
+*/
+void concatSparse(Eigen::SparseMatrix<double>& source,
+		std::vector<Triplet>& targetTriplets, int firstTargetRow) {
+	typedef Eigen::SparseMatrix<double>::InnerIterator Iterator;
+
+	// From Eigen tutorials...
+	for (int k = 0; k < source.outerSize(); k++) {
+		for (Iterator it(source, k); it; ++it) {
+			int sourceRow = it.row();
+			int column = it.col();
+			double value = it.value();
+			targetTriplets.push_back(Triplet(firstTargetRow + sourceRow, column, value));
+		}
+	}
+
+}
 
 void fd_grad(
-  const int nx,
-  const int ny,
-  const int nz,
-  const double h,
-  Eigen::SparseMatrix<double> & G)
-{
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+		const int nx,
+		const int ny,
+		const int nz,
+		const double h,
+		Eigen::SparseMatrix<double> & G) {
+
+	Eigen::SparseMatrix<double> Dx((nx - 1) * ny * nz, nx * ny * nz);
+	fd_partial_derivative(nx, ny, nz, h, 0, Dx);
+
+	Eigen::SparseMatrix<double> Dy(nx * (ny - 1) * nz, nx * ny * nz);
+	fd_partial_derivative(nx, ny, nz, h, 1, Dy);
+
+	Eigen::SparseMatrix<double> Dz(nx * ny * (nz - 1), nx * ny * nz);
+	fd_partial_derivative(nx, ny, nz, h, 2, Dz);
+
+	// Concatenate
+	std::vector<Triplet> triplets;
+	concatSparse(Dx, triplets, 0);
+	concatSparse(Dy, triplets, Dx.rows());
+	concatSparse(Dz, triplets, Dx.rows() + Dy.rows());	
+
+	G.resize(Dx.rows() + Dy.rows() + Dz.rows(), Dx.cols());
+	G.setFromTriplets(triplets.begin(), triplets.end());
+
 }

src/fd_interpolate.cpp

@@ -1,15 +1,61 @@
 #include "fd_interpolate.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)
-{
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+		const int nx,
+		const int ny,
+		const int nz,
+		const double h,
+		const Eigen::RowVector3d & corner,
+		const Eigen::MatrixXd & P,
+		Eigen::SparseMatrix<double> & W) {
+	W.resize(P.rows(), nx * ny * nz);
+
+	auto toGridIndex = [nx, ny, nz](int x, int y, int z) -> int {
+		return x + nx * y + nx * ny * z;
+	};
+
+	// As the Eigen guide says to, we build the matrix from triplets.
+	typedef Eigen::Triplet<double> Triplet;
+	std::vector<Triplet> tripletList;
+
+	int rowCount = P.rows();
+	for (int r = 0; r < rowCount; r++) {
+
+		// How many h's from the corner the point is.
+		double xFromStart = (P(r, 0) - corner(0)) / h;
+		double yFromStart = (P(r, 1) - corner(1)) / h;
+		double zFromStart = (P(r, 2) - corner(2)) / h;
+
+		// Indices of the lower corner of the grid cell.
+		int xGridIndex = floor(xFromStart);
+		int yGridIndex = floor(yFromStart);
+		int zGridIndex = floor(zFromStart);
+
+		// Where within the current grid cell the point is.
+		double x = xFromStart - xGridIndex;
+		double y = yFromStart - yGridIndex;
+		double z = zFromStart - zGridIndex;
+
+		// The weights computed by hand
+		double w0 = (1 - x) * (1 - z) * (1 - y);
+		double w1 = (1 - x) * z * (1 - y);
+		double w2 = (1 - x) * (1 - z) * y;
+		double w3 = (1 - x) * z * y;
+		double w4 = x * (1 - z) * (1 - y);
+		double w5 = x * z * (1 - y);
+		double w6 = x * (1 - z) * y;
+		double w7 = x * z * y;
+
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex, zGridIndex), w0));
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex, zGridIndex + 1), w1));
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex + 1, zGridIndex), w2));
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex, yGridIndex + 1, zGridIndex + 1), w3));
+
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex, zGridIndex), w4));
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex, zGridIndex + 1), w5));
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex + 1, zGridIndex), w6));
+		tripletList.push_back(Triplet(r, toGridIndex(xGridIndex + 1, yGridIndex + 1, zGridIndex + 1), w7));
+	}
+
+	W.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/fd_partial_derivative.cpp

@@ -1,14 +1,45 @@
 #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)
-{
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+		const int nx,
+		const int ny,
+		const int nz,
+		const double h,
+		const int dir,
+		Eigen::SparseMatrix<double> & D) {
+
+	// As the Eigen guide says to, we build the matrix from triplets.
+	typedef Eigen::Triplet<double> Triplet;
+	std::vector<Triplet> tripletList;
+
+	// Iterate over all cells. The limit depends on which direction we're
+	// going, so we can avoid going off the end.
+	int xLimit = dir == 0 ? nx - 1 : nx;
+	int yLimit = dir == 1 ? ny - 1 : ny;
+	int zLimit = dir == 2 ? nz - 1 : nz;
+	int r = 0;
+	for (int x = 0; x < xLimit; x++) {
+		for (int y = 0; y < yLimit; y++) {
+			for (int z = 0; z < zLimit; z++) {
+
+				// Co-efficient corresponding to one of the points in the
+				// finite difference. I divide by h to represent the division
+				// by dx.
+				int firstIndex = x + y * nx + z * nx * ny;
+				tripletList.push_back(Triplet(r, firstIndex, -1.0 / h));
+
+				// Second co-efficient. Again we have to divide by h.
+				int secondX = dir == 0 ? x + 1 : x;
+				int secondY = dir == 1 ? y + 1 : y;
+				int secondZ = dir == 2 ? z + 1 : z;
+				int secondIndex = secondX + secondY * nx + 
+						secondZ * nx * ny;
+				tripletList.push_back(Triplet(r, secondIndex, 1.0 / h));
+
+				r++;
+			}
+		}
+	}
+
+	D.setFromTriplets(tripletList.begin(), tripletList.end());
 }

src/poisson_surface_reconstruction.cpp

@@ -1,56 +1,106 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include "fd_interpolate.h"
+#include "fd_grad.h"
+#include<Eigen/IterativeLinearSolvers>
+
 
 void poisson_surface_reconstruction(
     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);
-
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
-
-  ////////////////////////////////////////////////////////////////////////////
-  // 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);
+	////////////////////////////////////////////////////////////////////////////
+	// 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);
+			}
+		}
+	}
+
+	/// Implementation Note:
+	/// Was not able to get it to correctly form the mesh. There must be a
+	/// bug somewhere, but despite long hours of searching, I have not been
+	/// able to solve it :(. With the exception of the bug, I am very
+	/// confident in my implementation.
+
+	// Weight matrices for each axis and the full weight matrix for the
+	// non-staggered grid.
+	Eigen::SparseMatrix<double> xWeightMatrix, yWeightMatrix, zWeightMatrix,
+			fullWeightMatrix;
+	fd_interpolate(nx - 1, ny, nz, h, corner + Eigen::RowVector3d(h / 2, 0, 0), P, xWeightMatrix);
+	fd_interpolate(nx, ny - 1, nz, h, corner + Eigen::RowVector3d(0, h /  2, 0), P, yWeightMatrix);
+	fd_interpolate(nx, ny, nz - 1, h, corner + Eigen::RowVector3d(0, 0, h / 2), P, zWeightMatrix);
+	fd_interpolate(nx, ny, nz, h, corner, P, fullWeightMatrix);
+
+	// Normals distributed onto grid.
+	Eigen::MatrixXd vX = xWeightMatrix.transpose() * N.col(0);
+	Eigen::MatrixXd vY = yWeightMatrix.transpose() * N.col(1);
+	Eigen::MatrixXd vZ = zWeightMatrix.transpose() * N.col(2);
+	Eigen::MatrixXd v(vX.rows() + vY.rows() + vZ.rows(), 1);
+	v << vX, vY, vZ;
+
+	// Gradient matrix via finite differences.
+	Eigen::SparseMatrix<double> gradient;
+	fd_grad(nx, ny, nz, h, gradient);
+
+	// Solve for the grid values. The standard form for least squares is 
+	// something like: Ax = b, so we have to compute the square gradient and
+	// apply the gradient transpose to v.
+	Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver;
+	Eigen::SparseMatrix<double> squareGradient = gradient.transpose() * gradient;
+	solver.compute(squareGradient);
+	Eigen::VectorXd g = solver.solve(gradient.transpose() * v);
+
+	// Find the iso value and shift g by it.
+	Eigen::MatrixXd sigmaProduct = fullWeightMatrix * g;
+	printf("Sigma product size: %d, %d\n", sigmaProduct.rows(), sigmaProduct.cols());
+	double sigma = sigmaProduct.sum() / n;
+	printf("Sigma: %f\n", sigma);
+	g.array() -= 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
+	////////////////////////////////////////////////////////////////////////////
+	igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F);
+
+	printf("Done\n");
 }