From 76e6d91cd03af05ac9d0b4a782ce5796522d7e45 Mon Sep 17 00:00:00 2001 From: Wenzheng Chen Date: Tue, 30 Oct 2018 23:15:23 -0400 Subject: [PATCH] new --- src/lscm.cpp | 60 +++++++++++++++++++++++++++---- src/tutte.cpp | 74 ++++++++++++++++++++++++++++++++++---- src/vector_area_matrix.cpp | 36 +++++++++++++++---- 3 files changed, 149 insertions(+), 21 deletions(-) diff --git a/src/lscm.cpp b/src/lscm.cpp index d6e1b6c..dca2ee4 100644 --- a/src/lscm.cpp +++ b/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 +#include +#include +#include +#include +#include +#include +#include + +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 L; + igl::cotmatrix(V, F, L); + + // A + Eigen::SparseMatrix A; + vector_area_matrix(F, A); + + // Q + Eigen::SparseMatrix LL; + LL = igl::repdiag(L, 2); + Eigen::SparseMatrix Q; + Q = L - A; + + // B + Eigen::SparseMatrix M; + igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M); + Eigen::SparseMatrix 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 svd(utu, + Eigen::ComputeFullU | Eigen::ComputeFullV); + U = U * svd.matrixU(); } diff --git a/src/tutte.cpp b/src/tutte.cpp index 6f6109b..cf2d304 100644 --- a/src/tutte.cpp +++ b/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 +#include +#include +#include +#include + +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 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& A, + const Eigen::MatrixBase & B, + const Eigen::MatrixBase & known, + const Eigen::MatrixBase & Y, + const Eigen::SparseMatrix& Aeq, + const Eigen::MatrixBase & Beq, + const bool pd, + Eigen::PlainObjectBase & Z + */ + + Eigen::SparseMatrix Aeq; + igl::min_quad_with_fixed_data 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); } diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp index 1d60ef8..483ced8 100644 --- a/src/vector_area_matrix.cpp +++ b/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& A) -{ - // Replace with your code - int V_size = F.maxCoeff()+1; - A.resize(V_size*2,V_size*2); +typedef Eigen::Triplet T; + +void vector_area_matrix(const Eigen::MatrixXi & F, + Eigen::SparseMatrix& 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 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 A_(V_size* 2, V_size*2); + A_.setFromTriplets(coef.begin(), coef.end()); + A = A_ + A_.transpose(); }