wenzheng chen hw1

src/fd_grad.cpp

@@ -1,13 +1,39 @@
 #include "fd_grad.h"
 
-void fd_grad(
-  const int nx,
-  const int ny,
-  const int nz,
-  const double h,
-  Eigen::SparseMatrix<double> & G)
-{
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+typedef Eigen::Triplet<double> T;
+
+void intervalue2(std::vector<T> &coef, int row, int aidx, int bidx) {
+	T tmp(row, aidx, -1);
+	T tmp2(row, bidx, 1);
+	coef.push_back(tmp);
+	coef.push_back(tmp2);
 }
+
+extern int gidx(double i, double j, double k, int nx, int ny, int nz);
+
+void fd_grad(const int nx, const int ny, const int nz, const double h,
+		Eigen::SparseMatrix<double> & G) {
+	////////////////////////////////////////////////////////////////////////////
+	// Add your code here
+
+	// first order
+	std::vector<T> coef;
+	for (int i = 0; i < nx - 1; i++)
+		for (int j = 0; j < ny - 1; j++)
+			for (int k = 0; k < nz - 1; k++) {
+
+				int st = gidx(i, j, k, nx, ny, nz);
+				int edx = gidx(i + 1, j, k, nx, ny, nz);
+				int edy = gidx(i, j + 1, k, nx, ny, nz);
+				int edz = gidx(i, j, k + 1, nx, ny, nz);
+
+				int row = st * 3;
+				intervalue2(coef, row, st, edx);
+				intervalue2(coef, row + 1, st, edy);
+				intervalue2(coef, row + 2, st, edz);
+			}
+
+	G.setFromTriplets(coef.begin(), coef.end());
+	////////////////////////////////////////////////////////////////////////////
+}
+

src/fd_interpolate.cpp

@@ -1,15 +1,86 @@
 #include "fd_interpolate.h"
+#include <assert.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
-  ////////////////////////////////////////////////////////////////////////////
+typedef Eigen::Triplet<double> T;
+
+void intervalue(std::vector<T> &coef, double aidx, double bidx, double value) {
+
+	for (int i = 0; i < 3; i++) {
+		T tmp(aidx * 3 + i, bidx * 3 + i, value);
+		coef.push_back(tmp);
+	}
+}
+
+double gidx(double i, double j, double k, int nx, int ny, int nz) {
+	double idx = i + nx * (j + k * ny);
+	assert(idx < (nx - 1) * (ny - 1) * (nz - 1));
+	return idx;
 }
+
+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
+
+	int pnum = P.rows();
+
+	// size of w should be pnum * 3 x (nx-1)(ny-1)(nz-1) * 3
+	std::vector<T> coef;
+	for (int p = 0; p < pnum; p++) {
+		double px = P(p, 0) - corner[0];
+		double py = P(p, 1) - corner[1];
+		double pz = P(p, 2) - corner[2];
+
+		double xind = px / h;
+		double yind = py / h;
+		double zind = pz / h;
+
+		// 8 equations
+		double xfloor = std::floor(px);
+		double yfloor = std::floor(py);
+		double zfloor = std::floor(pz);
+		double xceil = xfloor + 1;
+		double yceil = yfloor + 1;
+		double zceil = zfloor + 1;
+
+		// left, left, left
+		double idx = gidx(xfloor, yfloor, zfloor, nx, ny, nz);
+		double value = (xceil - px) * (yceil - py) * (zceil - pz);
+		intervalue(coef, p, idx, value);
+
+		idx = gidx(xfloor, yceil, zfloor, nx, ny, nz);
+		value = (xceil - px) * (py - yfloor) * (zceil - pz);
+		intervalue(coef, p, idx, value);
+
+		idx = gidx(xceil, yceil, zfloor, nx, ny, nz);
+		value = (px - xfloor) * (py - yfloor) * (zceil - pz);
+		intervalue(coef, p, idx, value);
+
+		idx = gidx(xceil, yfloor, zfloor, nx, ny, nz);
+		value = (px - xfloor) * (yceil - py) * (zceil - pz);
+		intervalue(coef, p, idx, value);
+
+		///////////////////////////////////////////////////////
+		idx = gidx(xfloor, yfloor, zceil, nx, ny, nz);
+		value = (xceil - px) * (yceil - px) * (pz - zfloor);
+		intervalue(coef, p, idx, value);
+
+		idx = gidx(xfloor, yceil, zceil, nx, ny, nz);
+		value = (xceil - px) * (px - yfloor) * (pz - zfloor);
+		intervalue(coef, p, idx, value);
+
+		idx = gidx(xceil, yceil, zceil, nx, ny, nz);
+		value = (px - xfloor) * (px - yfloor) * (pz - zfloor);
+		intervalue(coef, p, idx, value);
+
+		idx = gidx(xceil, yfloor, zceil, nx, ny, nz);
+		value = (px - xfloor) * (yceil - py) * (pz - zfloor);
+		intervalue(coef, p, idx, value);
+	}
+
+	W.setFromTriplets(coef.begin(), coef.end());
+
+	////////////////////////////////////////////////////////////////////////////
+}
+

src/fd_partial_derivative.cpp

@@ -1,14 +1,38 @@
 #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
-  ////////////////////////////////////////////////////////////////////////////
+typedef Eigen::Triplet<double> T;
+extern int gidx(double i, double j, double k, int nx, int ny, int nz);
+extern void intervalue2(std::vector<T> &coef, int row, int aidx, int bidx);
+
+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
+	// first order
+	std::vector<T> coef;
+	for (int i = 0; i < nx - 1; i++)
+		for (int j = 0; j < ny - 1; j++)
+			for (int k = 0; k < nz - 1; k++) {
+
+				int st = gidx(i, j, k, nx, ny, nz);
+				int edx = gidx(i + 1, j, k, nx, ny, nz);
+				int edy = gidx(i, j + 1, k, nx, ny, nz);
+				int edz = gidx(i, j, k + 1, nx, ny, nz);
+
+				int row = st;
+				switch (dir) {
+				case 0:
+					intervalue2(coef, row, st, edx);
+					break;
+				case 1:
+					intervalue2(coef, row + 1, st, edy);
+					break;
+				case 2:
+					intervalue2(coef, row + 2, st, edz);
+					break;
+				}
+			}
+
+	D.setFromTriplets(coef.begin(), coef.end());
+	////////////////////////////////////////////////////////////////////////////
 }

src/poisson_surface_reconstruction.cpp

@@ -2,55 +2,89 @@
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
 
-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);
+#include "fd_interpolate.h"
+#include "fd_partial_derivative.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
+	////////////////////////////////////////////////////////////////////////////
+	// first, we calculate spaese mtx M, which M_p3_xyz3 * X_(x-1)(y-1)(z-1)3_1 = P_p3_1
+	int pnum = n;
+	int gnum = (nx - 1) * (ny - 1) * (nz - 1);
+	Eigen::SparseMatrix<double> W(pnum * 3, gnum * 3);
+	fd_interpolate(nx, ny, nz, h, corner, P, W);
+
+	// next, we calculate sparse mtx G, which G_(x-1)(y-1)(z-1)3_xyz * X_xyz_1 = D_(x-1)(y-1)(z-1)_1
+	Eigen::SparseMatrix<double> G((nx - 1) * (ny - 1) * (nz - 1) * 3,
+			nx * ny * nz * 1);
+	fd_grad(nx, ny, nz, h, G);
+
+	// next, we combine G with M
+	Eigen::SparseMatrix<double> A = W * G;
+
+	// last, we solve Ax = B
+	Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
+	solver.compute(A);
+
+	Eigen::VectorXd b(n * 3);
+	for (int i = 0; i < pnum; i++) {
+		for (int j = 0; j < 3; j++) {
+			b(i * 3 + j, 0) = N(i, j);
+		}
+	}
+	g = solver.solve(b);
+
+	double sigma = (W * g).mean();
+	for (int i = 0; i < g.cols(); i++)
+		g(i) = g(i) - 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);
 }