From 6a230d3e7e4f081ab60c16a9ab8f91005da5d54e Mon Sep 17 00:00:00 2001 From: Adrian Date: Mon, 29 Oct 2018 00:38:33 -0400 Subject: [PATCH 1/3] tutte coordinates ok, weird orientation flips happening --- src/tutte.cpp | 69 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 68 insertions(+), 1 deletion(-) diff --git a/src/tutte.cpp b/src/tutte.cpp index 6f6109b..725228f 100644 --- a/src/tutte.cpp +++ b/src/tutte.cpp @@ -1,4 +1,9 @@ #include "tutte.h" +#include "igl/boundary_loop.h" +#include "igl/map_vertices_to_circle.h" +#include "igl/edges.h" +#include "igl/min_quad_with_fixed.h" +#include void tutte( const Eigen::MatrixXd & V, @@ -6,6 +11,68 @@ void tutte( Eigen::MatrixXd & U) { // Replace with your code - U = V.leftCols(2); + + // Determine all boundary loop on the mesh + std::vector> bdL; + igl::boundary_loop(F, bdL); + + // Identify the largest boundary loop + int max_idx = 0; + for (int i = 0; i < bdL.size(); i++) { + if (bdL[i].size() >= bdL[max_idx].size()) { + max_idx = i; + } + } + Eigen::VectorXi bdry(bdL[max_idx].size()); + for (int i = 0; i < bdL[max_idx].size(); i++) { + bdry[i] = bdL[max_idx][i]; + } + std::cout << "CHANGED BOUNDARY MAP TO AN EIGEN::VECTOR" << std::endl; + Eigen::MatrixXd bdry_pos; + igl::map_vertices_to_circle(V, bdry, bdry_pos); + std::cout << "OK figured out what the boundary is" << std::endl; + + // Now apply the weighted Tutte map + Eigen::MatrixXi E; + igl::edges(F, E); + std::cout << "OK figured out what the edges are too " << std::endl; + + Eigen::SparseMatrix L; + typedef Eigen::Triplet T; + std::vector triples; + triples.reserve(4 * E.rows()); + for (int i = 0; i < E.rows(); i++) { + double weight = 1.0 / ((V.row(E(i,0)) - V.row(E(i,1))).norm()); + triples.push_back(T(E(i,0), E(i,1), weight)); + triples.push_back(T(E(i,1), E(i,0), weight)); + triples.push_back(T(E(i,0), E(i,0),-weight)); + triples.push_back(T(E(i,1), E(i,1),-weight)); + } + int n_v = F.maxCoeff() + 1; + L.resize(n_v, n_v); + L.setFromTriplets(triples.begin(), triples.end()); + std::cout << "figured out the Laplacian" << std::endl; + + // Now solve the weighted optimization problem subject to boundary constraints U_x^T L U_x + U_y^T L U_y + igl::min_quad_with_fixed_data data; + + // No linear constraints or linear term in optimization objective; + Eigen::SparseMatrix Z; + Z.setZero(); + Eigen::VectorXd B; + B.setZero(); + Eigen::VectorXd z(n_v); + z.setZero(); + + igl::min_quad_with_fixed_precompute(L, bdry, Z, false, data); // no linear equality constraints + std::cout << "factorized data out laplacian" << std::endl; + U.resize(n_v, 2); + for (int i = 0; i <= 1; i++) { + Eigen::VectorXd cur_col; + igl::min_quad_with_fixed_solve(data, z, bdry_pos.col(i), B, cur_col); // no linear term in objective or linear inequality contraints + std::cout << "I SOLVED IT " << i << std::endl; + U.col(i) = cur_col; + } + std::cout << "figured out full paramterization" << std::endl; } From ea795acedf4bb8af223e5115809f5d0556f82677 Mon Sep 17 00:00:00 2001 From: Adrian Date: Mon, 29 Oct 2018 22:53:55 -0400 Subject: [PATCH 2/3] submit --- src/lscm.cpp | 53 +++++++++++++++++++++++++++++++++++--- src/tutte.cpp | 33 ++++-------------------- src/vector_area_matrix.cpp | 26 ++++++++++++++++++- 3 files changed, 80 insertions(+), 32 deletions(-) diff --git a/src/lscm.cpp b/src/lscm.cpp index d6e1b6c..36c4c03 100644 --- a/src/lscm.cpp +++ b/src/lscm.cpp @@ -1,10 +1,57 @@ #include "lscm.h" +#include "vector_area_matrix.h" +#include "igl/repdiag.h" +#include "igl/cotmatrix.h" +#include "igl/massmatrix.h" +#include "igl/eigs.h" +#include void lscm( const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::MatrixXd & U) -{ - // Replace with your code - U = V.leftCols(2); +{ + // Construct matrix Q = (L 0; 0 L) - A where L is the laplacian and A is the vector area + Eigen::SparseMatrix L; + igl::cotmatrix(V, F, L); + Eigen::SparseMatrix L_d; + igl::repdiag(L, 2, L_d); + Eigen::SparseMatrix A; + vector_area_matrix(F, A); + Eigen::SparseMatrix Q = L_d - A; + + // Construct the matrix B = (M 0; 0 M) where M is the mass matrix + Eigen::SparseMatrix M; + igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_BARYCENTRIC, M); + Eigen::SparseMatrix B; + igl::repdiag(M, 2, B); + + // Find parameterization coordinates by a eigenvector of Qv = lambda * B * v with v non-zero eigenvalue + Eigen::VectorXd eigen_vals; + Eigen::MatrixXd eigen_vecs; + igl::eigs(Q, B, 5, igl::EigsType::EIGS_TYPE_SM, eigen_vecs, eigen_vals); + + int idx = 1; // the second eigenvector calculated should be the Fielder vector + double tol = 1e-4; // Pick an eigenvector sufficiently away from 0 + // hope that an eigen_vec with non-zero eigen_val has been found among the calculated evals + // otherwise just take the largest evec with the largest evec. + while (abs(eigen_vals[idx]) <= tol) { + idx++; + } + if (idx >= 5) { + idx = 4; + } + + // Reshape eigenvector into U + int n_v = F.maxCoeff() + 1; + U.resize(n_v, 2); + U.col(0) = eigen_vecs.col(idx).head(n_v); + U.col(1) = eigen_vecs.col(idx).tail(n_v); + + // Do SVD on coordinates U to find a rotation on the coordinates U + Eigen::Matrix2d covar = U.transpose() * U; + Eigen::JacobiSVD svd(covar, Eigen::ComputeFullU); + for (int i = 0; i < U.rows(); i++) { + U.row(i) = U.row(i) * svd.matrixU().transpose(); + } } diff --git a/src/tutte.cpp b/src/tutte.cpp index 725228f..371beb5 100644 --- a/src/tutte.cpp +++ b/src/tutte.cpp @@ -3,40 +3,21 @@ #include "igl/map_vertices_to_circle.h" #include "igl/edges.h" #include "igl/min_quad_with_fixed.h" -#include void tutte( const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::MatrixXd & U) { - // Replace with your code - - // Determine all boundary loop on the mesh - std::vector> bdL; - igl::boundary_loop(F, bdL); - - // Identify the largest boundary loop - int max_idx = 0; - for (int i = 0; i < bdL.size(); i++) { - if (bdL[i].size() >= bdL[max_idx].size()) { - max_idx = i; - } - } - Eigen::VectorXi bdry(bdL[max_idx].size()); - for (int i = 0; i < bdL[max_idx].size(); i++) { - bdry[i] = bdL[max_idx][i]; - } - std::cout << "CHANGED BOUNDARY MAP TO AN EIGEN::VECTOR" << std::endl; + // Determine the longest boundary loop and map its vertices onto the unit circle + Eigen::VectorXi bdry; + igl::boundary_loop(F, bdry); Eigen::MatrixXd bdry_pos; igl::map_vertices_to_circle(V, bdry, bdry_pos); - std::cout << "OK figured out what the boundary is" << std::endl; - // Now apply the weighted Tutte map + // Calculate the weighted Laplacian matrix Eigen::MatrixXi E; igl::edges(F, E); - std::cout << "OK figured out what the edges are too " << std::endl; - Eigen::SparseMatrix L; typedef Eigen::Triplet T; std::vector triples; @@ -51,9 +32,8 @@ void tutte( int n_v = F.maxCoeff() + 1; L.resize(n_v, n_v); L.setFromTriplets(triples.begin(), triples.end()); - std::cout << "figured out the Laplacian" << std::endl; - // Now solve the weighted optimization problem subject to boundary constraints U_x^T L U_x + U_y^T L U_y + // Now solve the weighted optimization problem U_x^T L U_x + U_y^T L U_y subject to boundary constraints calculated earlier igl::min_quad_with_fixed_data data; // No linear constraints or linear term in optimization objective; @@ -65,14 +45,11 @@ void tutte( z.setZero(); igl::min_quad_with_fixed_precompute(L, bdry, Z, false, data); // no linear equality constraints - std::cout << "factorized data out laplacian" << std::endl; U.resize(n_v, 2); for (int i = 0; i <= 1; i++) { Eigen::VectorXd cur_col; igl::min_quad_with_fixed_solve(data, z, bdry_pos.col(i), B, cur_col); // no linear term in objective or linear inequality contraints - std::cout << "I SOLVED IT " << i << std::endl; U.col(i) = cur_col; } - std::cout << "figured out full paramterization" << std::endl; } diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp index 1d60ef8..ffa58da 100644 --- a/src/vector_area_matrix.cpp +++ b/src/vector_area_matrix.cpp @@ -1,11 +1,35 @@ #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); + + // Determine the longest boundary loop + Eigen::VectorXi bdry; + igl::boundary_loop(F, bdry); + + // Each edge (i, j) in the loop contributes 1/2 (x_i y_j - x_j y_i) to the vector area + typedef Eigen::Triplet T; + std::vector triples; + triples.reserve(bdry.size() * 4); + for (int i = 0; i < bdry.size(); i++) { + if (i == bdry.size() - 1) { + triples.push_back(T(bdry[i], bdry[0] + V_size, 0.25)); + triples.push_back(T(bdry[0] + V_size, bdry[i], 0.25)); + triples.push_back(T(bdry[0], bdry[i] + V_size, -0.25)); + triples.push_back(T(bdry[i] + V_size, bdry[0], -0.25)); + } + else { + triples.push_back(T(bdry[i], bdry[i+1] + V_size, 0.25)); + triples.push_back(T(bdry[i+1] + V_size, bdry[i], 0.25)); + triples.push_back(T(bdry[i+1], bdry[i] + V_size, -0.25)); + triples.push_back(T(bdry[i] + V_size, bdry[i+1], -0.25)); + } + } + A.setFromTriplets(triples.begin(), triples.end()); } From 66eb87e2893579f8cc1e545286a2ee851800333f Mon Sep 17 00:00:00 2001 From: Adrian Date: Tue, 30 Oct 2018 17:30:54 -0400 Subject: [PATCH 3/3] fix lscm method weirdness --- src/lscm.cpp | 6 +++--- src/vector_area_matrix.cpp | 29 ++++++++++++++++------------- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/lscm.cpp b/src/lscm.cpp index 36c4c03..1d11d9c 100644 --- a/src/lscm.cpp +++ b/src/lscm.cpp @@ -18,7 +18,7 @@ void lscm( igl::repdiag(L, 2, L_d); Eigen::SparseMatrix A; vector_area_matrix(F, A); - Eigen::SparseMatrix Q = L_d - A; + Eigen::SparseMatrix Q = L_d - 2.0 * A; // Construct the matrix B = (M 0; 0 M) where M is the mass matrix Eigen::SparseMatrix M; @@ -32,7 +32,7 @@ void lscm( igl::eigs(Q, B, 5, igl::EigsType::EIGS_TYPE_SM, eigen_vecs, eigen_vals); int idx = 1; // the second eigenvector calculated should be the Fielder vector - double tol = 1e-4; // Pick an eigenvector sufficiently away from 0 + double tol = 1e-8; // Pick an eigenvector sufficiently away from 0 // hope that an eigen_vec with non-zero eigen_val has been found among the calculated evals // otherwise just take the largest evec with the largest evec. while (abs(eigen_vals[idx]) <= tol) { @@ -52,6 +52,6 @@ void lscm( Eigen::Matrix2d covar = U.transpose() * U; Eigen::JacobiSVD svd(covar, Eigen::ComputeFullU); for (int i = 0; i < U.rows(); i++) { - U.row(i) = U.row(i) * svd.matrixU().transpose(); + U.row(i) = U.row(i) * svd.matrixU(); } } diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp index ffa58da..2195b99 100644 --- a/src/vector_area_matrix.cpp +++ b/src/vector_area_matrix.cpp @@ -6,30 +6,33 @@ void vector_area_matrix( Eigen::SparseMatrix& A) { int V_size = F.maxCoeff()+1; - A.resize(V_size*2,V_size*2); + Eigen::SparseMatrix A_u; + A_u.resize(V_size*2,V_size*2); // Determine the longest boundary loop - Eigen::VectorXi bdry; - igl::boundary_loop(F, bdry); + std::vector> bdry_loops; + igl::boundary_loop(F, bdry_loops); // Each edge (i, j) in the loop contributes 1/2 (x_i y_j - x_j y_i) to the vector area typedef Eigen::Triplet T; std::vector triples; - triples.reserve(bdry.size() * 4); + triples.reserve(bdry_loops[0].size() * 4); + for (int l = 0; l < bdry_loops.size() ; l++) { + std::vector bdry = bdry_loops[l]; for (int i = 0; i < bdry.size(); i++) { if (i == bdry.size() - 1) { - triples.push_back(T(bdry[i], bdry[0] + V_size, 0.25)); - triples.push_back(T(bdry[0] + V_size, bdry[i], 0.25)); - triples.push_back(T(bdry[0], bdry[i] + V_size, -0.25)); - triples.push_back(T(bdry[i] + V_size, bdry[0], -0.25)); + triples.push_back(T(bdry[i], bdry[0] + V_size, 0.5)); + triples.push_back(T(bdry[i] + V_size, bdry[0], -0.5)); } else { - triples.push_back(T(bdry[i], bdry[i+1] + V_size, 0.25)); - triples.push_back(T(bdry[i+1] + V_size, bdry[i], 0.25)); - triples.push_back(T(bdry[i+1], bdry[i] + V_size, -0.25)); - triples.push_back(T(bdry[i] + V_size, bdry[i+1], -0.25)); + triples.push_back(T(bdry[i], bdry[i+1] + V_size, 0.5)); + triples.push_back(T(bdry[i] + V_size, bdry[i+1], -0.5)); } } - A.setFromTriplets(triples.begin(), triples.end()); +} + A_u.setFromTriplets(triples.begin(), triples.end()); + Eigen::SparseMatrix B; + B = A_u.transpose(); + A = (A_u + B) / 2; }