wenzheng hw 5

src/lscm.cpp

@@ -1,10 +1,56 @@
 #include "lscm.h"
 
-void lscm(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXd & U)
-{
-  // Replace with your code
-  U = V.leftCols(2);
+#include "vector_area_matrix.h"
+
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/repdiag.h>
+#include <igl/massmatrix.h>
+#include <igl/eigs.h>
+
+void lscm(const Eigen::MatrixXd & V, const Eigen::MatrixXi & F,
+		Eigen::MatrixXd & U) {
+	// Replace with your code
+	//  U = V.leftCols(2);
+
+	// L
+	// calculate L, here, w is 1
+	Eigen::SparseMatrix<double> L;
+	igl::cotmatrix(V, F, L);
+
+	// A
+	Eigen::SparseMatrix<double> A;
+	vector_area_matrix(F, A);
+
+	// Q
+	Eigen::SparseMatrix<double> LL;
+	LL = igl::repdiag(L, 2);
+	Eigen::SparseMatrix<double> Q;
+	Q = L - A;
+
+	// B
+	Eigen::SparseMatrix<double> M;
+	igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+	Eigen::SparseMatrix<double> B;
+	B = igl::repdiag(M, 2);
+
+	// equation
+	Eigen::MatrixXd Vecs;
+	Eigen::VectorXd lambdas;
+	igl::eigs(Q, B, 3, igl::EIGS_TYPE_SM, Vecs, lambdas);
+
+	int vnum = V.rows();
+	U.resize(vnum, 2);
+	U.col(0) = Vecs.col(2).head(vnum);
+	U.col(1) = Vecs.col(2).tail(vnum);
+
+	// final rotation
+	Eigen::MatrixXd utu = U.transpose() * U;
+	// svd
+	Eigen::JacobiSVD<Eigen::MatrixXd> svd(utu,
+			Eigen::ComputeFullU | Eigen::ComputeFullV);
+	U = U * svd.matrixU();
 }

src/tutte.cpp

@@ -1,11 +1,71 @@
 #include "tutte.h"
 
-void tutte(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXd & U)
-{
-  // Replace with your code
-  U = V.leftCols(2);
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/min_quad_with_fixed.h>
+
+void tutte(const Eigen::MatrixXd & V, const Eigen::MatrixXi & F,
+		Eigen::MatrixXd & U) {
+	// Replace with your code
+	// U = V.leftCols(2);
+
+	int vnum = V.rows();
+	int fnum = F.rows()
+
+	// calculate L, here, w is 1
+	Eigen::SparseMatrix<double> L;
+	igl::cotmatrix(V, F, L);
+
+	// map the boundary into circle
+	// we have some known and some unknown
+	Eigen::VectorXd bounds;
+	igl::boundary_loop(F, bounds);
+
+	// these points will be maped into circles
+	Eigen::MatrixXd UVknown;
+	igl::map_vertices_to_circle(V, bounds, UVknown);
+
+	/*
+	 *   // MIN_QUAD_WITH_FIXED Minimize a quadratic energy of the form
+	 //
+	 // trace( 0.5*Z'*A*Z + Z'*B + constant )
+	 //
+	 // subject to
+	 //
+	 //   Z(known,:) = Y, and
+	 //   Aeq*Z = Beq
+	 //   T  should be a eigen matrix primitive type like int or double
+	 // Inputs:
+	 //   A  n by n matrix of quadratic coefficients
+	 //   known list of indices to known rows in Z
+	 //   Y  list of fixed values corresponding to known rows in Z
+	 //   Aeq  m by n list of linear equality constraint coefficients
+	 //   pd flag specifying whether A(unknown,unknown) is positive definite
+	 // Outputs:
+	 //   data  factorization struct with all necessary information to solve
+	 //     using min_quad_with_fixed_solve
+	 // Returns true on success, false on error
+	 */
+	// based on these knows, we solve unknows by minimize trace(UtLU)
+	/*
+	 *   const Eigen::SparseMatrix<T>& A,
+	 const Eigen::MatrixBase<DerivedB> & B,
+	 const Eigen::MatrixBase<Derivedknown> & known,
+	 const Eigen::MatrixBase<DerivedY> & Y,
+	 const Eigen::SparseMatrix<T>& Aeq,
+	 const Eigen::MatrixBase<DerivedBeq> & Beq,
+	 const bool pd,
+	 Eigen::PlainObjectBase<DerivedZ> & Z
+	 */
+
+	Eigen::SparseMatrix<double> Aeq;
+	igl::min_quad_with_fixed_data<double> data;
+	min_quad_with_fixed_precompute(L, bounds, Aeq, false, data);
+
+	Eigen::VectorXd B = Eigen::VectorXd::Zero(vnum);
+	Eigen::VectorXd Beq = Eigen::VectorXd::Zero(vnum);
+	min_quad_with_fixed_solve(data, B, UVknown, Beq, U);
 }
 

src/vector_area_matrix.cpp

@@ -1,11 +1,33 @@
 #include "vector_area_matrix.h"
+#include "igl/boundary_loop.h"
 
-void vector_area_matrix(
-  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);
+typedef Eigen::Triplet<double> T;
+
+void vector_area_matrix(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);
+
+	Eigen::VectorXd bounds;
+	igl::boundary_loop(F, bounds);
+	int bnum = bounds.size();
+
+	std::vector<T> coef;
+	for (int i = 0; i < bnum; i++) {
+		int b1 = bounds[i];
+		for (int j = 0; j < bnum; j++) {
+			int b2 = bounds[j];
+			T tmp(b1, b2 + V_size, 0.5);
+			coef.push_back(tmp);
+
+			T tmp2(b1 + V_size, b2, -0.5);
+			coef.push_back(tmp2);
+		}
+	}
+
+	Eigen::SparseMatrix<double> A_(V_size* 2, V_size*2);
+	A_.setFromTriplets(coef.begin(), coef.end());
+	A = A_ + A_.transpose();
 }