Zhiyuan Cheng's HW05

src/lscm.cpp

@@ -1,10 +1,35 @@
 #include "lscm.h"
+#include "vector_area_matrix.h"
+#include <igl/cotmatrix.h>
+#include <igl/repdiag.h>
+#include <igl/massmatrix.h>
+#include <igl/eigs.h>
+#include <Eigen/SVD>
 
 void lscm(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXd & U)
+	const Eigen::MatrixXd & V,
+	const Eigen::MatrixXi & F,
+	Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
+	// construct Q and B
+	Eigen::SparseMatrix<double> L, L2, A, Q, M, B;
+	igl::cotmatrix(V, F, L);
+	igl::repdiag(L, 2, L2);
+	vector_area_matrix(F, A);
+	Q.resize(A.rows(), A.rows()); // A is 2#V x 2#V
+	Q = L2 - 2.0 * A;
+	igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_DEFAULT, M);
+	igl::repdiag(M, 2, B);
+
+	// solve generalized Eigen value problem
+	Eigen::MatrixXd eigenVec;
+	Eigen::VectorXd eigenVal;
+	igl::eigs(Q, B, 3, igl::EigsType::EIGS_TYPE_SM, eigenVec, eigenVal);
+
+	U.resize(V.rows(), 2);
+	U.col(0) = eigenVec.col(2).head(V.rows());
+	U.col(1) = eigenVec.col(2).tail(V.rows());
+	// rotate
+	Eigen::JacobiSVD<Eigen::MatrixXd> svd(U.transpose() * U, Eigen::ComputeFullU | Eigen::ComputeFullV);
+	U = U * svd.matrixU();
 }

src/tutte.cpp

@@ -1,11 +1,47 @@
 #include "tutte.h"
+#include <igl/cotmatrix.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/boundary_loop.h>
+#include <igl/min_quad_with_fixed.h>
+
 
 void tutte(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXd & U)
+	const Eigen::MatrixXd & V,
+	const Eigen::MatrixXi & F,
+	Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
-}
+	Eigen::VectorXi boundary;
+	igl::boundary_loop(F, boundary);
+
+	Eigen::MatrixXd UV;
+	igl::map_vertices_to_circle(V, boundary, UV);
 
+	Eigen::SparseMatrix<double> L;
+	// igl::cotmatrix(V, F, L);
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	for (int f = 0; f < F.rows(); f++) {
+		for (int i = 0; i < 3; i++) {
+			int j = (i + 1) % 3;
+			int k = (i + 2) % 3;
+			double norm_ij = (V.row(F(f, i)) - V.row(F(f, j))).norm();
+			double norm_ik = (V.row(F(f, i)) - V.row(F(f, k))).norm();
+
+			// case i != j
+			tripletList.push_back(T(F(f, i), F(f, j), 1.0/norm_ij));
+			tripletList.push_back(T(F(f, i), F(f, k), 1.0/norm_ik));
+			// case i == j
+			tripletList.push_back(T(F(f, i), F(f, i), -1.0/norm_ij));
+			tripletList.push_back(T(F(f, i), F(f, i), -1.0/norm_ik));
+		}
+	}
+	L.resize(V.rows(), V.rows());
+	L.setFromTriplets(tripletList.begin(), tripletList.end());
+
+	Eigen::SparseMatrix<double> Aeq;
+	Eigen::VectorXd B = Eigen::VectorXd::Zero(V.rows());
+	igl::min_quad_with_fixed(L, B, boundary, UV, Aeq, B, false, U);
+
+	// igl::min_quad_with_fixed_precompute(L, boundary, Aeq, false, data);
+	// igl::min_quad_with_fixed_solve(data, B, UV, Beq, U);
+}

src/vector_area_matrix.cpp

@@ -1,11 +1,28 @@
 #include "vector_area_matrix.h"
+#include <igl/boundary_loop.h>
 
 void vector_area_matrix(
-  const Eigen::MatrixXi & F,
-  Eigen::SparseMatrix<double>& A)
+	const Eigen::MatrixXi & F,
+	Eigen::SparseMatrix<double>& A)
 {
-  // Replace with your code
-  int V_size = F.maxCoeff()+1;
-  A.resize(V_size*2,V_size*2);
+	int V_size = F.maxCoeff() + 1;
+
+	std::vector<std::vector<int>> boundary;
+	igl::boundary_loop(F, boundary);
+
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	for (int b = 0; b < boundary.size(); b++) {
+		std::vector<int> bound = boundary[b];
+		for (int i = 0; i < bound.size(); i++) {
+			int j = (i + 1) % bound.size();
+			tripletList.push_back(T(bound[i], bound[j] + V_size, 0.25));
+			tripletList.push_back(T(bound[i] + V_size, bound[j], -0.25));
+			tripletList.push_back(T(bound[j], bound[i] + V_size, -0.25));
+			tripletList.push_back(T(bound[j] + V_size, bound[i], 0.25));
+		}
+	}
+	A.resize(V_size * 2, V_size * 2);
+	A.setFromTriplets(tripletList.begin(), tripletList.end());
 }