From 1b43febb83e1590880d37e59e97ee623424345da Mon Sep 17 00:00:00 2001 From: Zhiyuan Date: Mon, 29 Oct 2018 17:42:32 -0400 Subject: [PATCH 1/2] submit --- src/lscm.cpp | 34 +++++++++++++++++++++++++++++----- src/tutte.cpp | 29 +++++++++++++++++++++++------ src/vector_area_matrix.cpp | 28 +++++++++++++++++++++++----- 3 files changed, 75 insertions(+), 16 deletions(-) diff --git a/src/lscm.cpp b/src/lscm.cpp index d6e1b6c..eedd96e 100644 --- a/src/lscm.cpp +++ b/src/lscm.cpp @@ -1,10 +1,34 @@ #include "lscm.h" +#include "vector_area_matrix.h" +#include +#include +#include +#include +#include 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 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 - 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 << eigenVec.col(2); + // rotate + Eigen::JacobiSVD svd(U.transpose() * U, Eigen::ComputeFullU | Eigen::ComputeFullV); + U = U * svd.matrixU(); } diff --git a/src/tutte.cpp b/src/tutte.cpp index 6f6109b..274a7e7 100644 --- a/src/tutte.cpp +++ b/src/tutte.cpp @@ -1,11 +1,28 @@ #include "tutte.h" +#include +#include +#include +#include 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 L; + igl::cotmatrix(V, F, L); + + Eigen::SparseMatrix Aeq; + Eigen::MatrixXd B, Beq; + B.resize(V.rows(), 2); + igl::min_quad_with_fixed(L, B, boundary, UV, Aeq, Beq, false, U); + + // igl::min_quad_with_fixed_precompute(L, boundary, Aeq, false, data); + // igl::min_quad_with_fixed_solve(data, B, UV, Beq, U); +} diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp index 1d60ef8..fb32a69 100644 --- a/src/vector_area_matrix.cpp +++ b/src/vector_area_matrix.cpp @@ -1,11 +1,29 @@ #include "vector_area_matrix.h" +#include void vector_area_matrix( - const Eigen::MatrixXi & F, - Eigen::SparseMatrix& A) + 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); + int V_size = F.maxCoeff() + 1; + + std::vector> boundary; + igl::boundary_loop(F, boundary); + + typedef Eigen::Triplet T; + std::vector tripletList; + tripletList.reserve(V_size * V_size * 4); + for (int b = 0; b < boundary.size(); b++) { + std::vector 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()); } From 40971d6d494a3f44dbc670b3c510f3ca9a1a5f2c Mon Sep 17 00:00:00 2001 From: Zhiyuan Date: Tue, 30 Oct 2018 21:44:54 -0400 Subject: [PATCH 2/2] fix L for tutte, fix eigs not converge --- src/lscm.cpp | 7 ++++--- src/tutte.cpp | 27 +++++++++++++++++++++++---- src/vector_area_matrix.cpp | 5 ++--- 3 files changed, 29 insertions(+), 10 deletions(-) diff --git a/src/lscm.cpp b/src/lscm.cpp index eedd96e..7d24785 100644 --- a/src/lscm.cpp +++ b/src/lscm.cpp @@ -17,7 +17,7 @@ void lscm( 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 - A; + Q = L2 - 2.0 * A; igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_DEFAULT, M); igl::repdiag(M, 2, B); @@ -25,9 +25,10 @@ void lscm( Eigen::MatrixXd eigenVec; Eigen::VectorXd eigenVal; igl::eigs(Q, B, 3, igl::EigsType::EIGS_TYPE_SM, eigenVec, eigenVal); - + U.resize(V.rows(), 2); - U << eigenVec.col(2); + U.col(0) = eigenVec.col(2).head(V.rows()); + U.col(1) = eigenVec.col(2).tail(V.rows()); // rotate Eigen::JacobiSVD svd(U.transpose() * U, Eigen::ComputeFullU | Eigen::ComputeFullV); U = U * svd.matrixU(); diff --git a/src/tutte.cpp b/src/tutte.cpp index 274a7e7..06425e2 100644 --- a/src/tutte.cpp +++ b/src/tutte.cpp @@ -4,6 +4,7 @@ #include #include + void tutte( const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, @@ -16,12 +17,30 @@ void tutte( igl::map_vertices_to_circle(V, boundary, UV); Eigen::SparseMatrix L; - igl::cotmatrix(V, F, L); + // igl::cotmatrix(V, F, L); + typedef Eigen::Triplet T; + std::vector 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 Aeq; - Eigen::MatrixXd B, Beq; - B.resize(V.rows(), 2); - igl::min_quad_with_fixed(L, B, boundary, UV, Aeq, Beq, false, U); + 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); diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp index fb32a69..277882e 100644 --- a/src/vector_area_matrix.cpp +++ b/src/vector_area_matrix.cpp @@ -12,15 +12,14 @@ void vector_area_matrix( typedef Eigen::Triplet T; std::vector tripletList; - tripletList.reserve(V_size * V_size * 4); for (int b = 0; b < boundary.size(); b++) { std::vector 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)); + 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);