From a50ea27a371334d18c1ca40a77e953b510c9af2a Mon Sep 17 00:00:00 2001 From: Wenzheng Chen Date: Sat, 13 Oct 2018 20:25:39 -0400 Subject: [PATCH] jw3 --- main.cpp | 12 ++--- src/cotmatrix.cpp | 106 ++++++++++++++++++++++++++++++++++++++++++--- src/massmatrix.cpp | 50 ++++++++++++++++++--- src/smooth.cpp | 41 ++++++++++++++---- 4 files changed, 182 insertions(+), 27 deletions(-) diff --git a/main.cpp b/main.cpp index c1d5fed..acca1a7 100644 --- a/main.cpp +++ b/main.cpp @@ -2,7 +2,7 @@ #include #include #include -#include +#include #include #include #include @@ -31,7 +31,7 @@ int main(int argc, char *argv[]) } } - igl::viewer::Viewer viewer; + igl::opengl::glfw::Viewer viewer; std::cout< T; -void cotmatrix( - const Eigen::MatrixXd & l, - const Eigen::MatrixXi & F, - Eigen::SparseMatrix & L) -{ - // Add your code here +void setV2F(Eigen::MatrixXi & V2F1, Eigen::MatrixXi & V2F2, int i, int j, + int k) { + + if (V2F1(i, j) == -1) { + V2F1(i, j) = k; + V2F1(j, i) = k; + } else { + V2F2(i, j) = k; + V2F2(j, i) = k; + } + return; +} + +void cotmatrix(const Eigen::MatrixXd & l, const Eigen::MatrixXi & F, + Eigen::SparseMatrix & L) { + // Add your code here + + int vnum = F.maxCoeff() + 1; + int fnum = F.rows(); + + // calculate area of faces + Eigen::VectorXd FA(fnum); + for (int i = 0; i < fnum; i++) { + int v1 = F(i, 0); + int v2 = F(i, 1); + int v3 = F(i, 2); + double e1len = l(v1, v2); + double e2len = l(v1, v3); + double e3len = l(v2, v3); + double s = e1len + e2len + e3len; + s = s / 2; + double A2 = s * (s - e1len) * (s - e2len) * (s - e3len); + double A = sqrt(A2); + FA[i] = A; + } + + // build point2face connection + // 2 point i, j will link to 2 faces + Eigen::MatrixXi V2F1, V2F2; + V2F1 = -Eigen::MatrixXi::Ones(vnum, vnum); + V2F2 = -Eigen::MatrixXi::Ones(vnum, vnum); + for (int i = 0; i < fnum; i++) { + int v1 = F(i, 0); + int v2 = F(i, 1); + int v3 = F(i, 2); + setV2F(V2F1, V2F2, v1, v2, i); + } + + std::vector coef; + for (int i = 0; i < vnum; i++) { + // diag + double ii = 0; + for (int j = 0; j < vnum; j++) { + + if (i == j) + continue; + + // is i and j connected? + int k1 = V2F1(i, j); + int k2 = V2F2(i, j); + // cot c = (a^2 + b^2 - c^2) / 4A + + int k = k1; + if (k != -1) { + int v1 = F(k, 0); + int v2 = F(k, 1); + int v3 = F(k, 2); + double e1len = l(v1, v2); + double e2len = l(v1, v3); + double e3len = l(v2, v3); + double A = FA[k]; + double cot = e1len * e1len + e2len * e2len - e3len * e3len; + cot = cot / (4 * A); + T tmp(i, j, cot); + coef.push_back(tmp); + ii += cot; + } + + k = k2; + if (k != -1) { + int v1 = F(k, 0); + int v2 = F(k, 1); + int v3 = F(k, 2); + double e1len = l(v1, v2); + double e2len = l(v1, v3); + double e3len = l(v2, v3); + double A = FA[k]; + double cot = e1len * e1len + e2len * e2len - e3len * e3len; + cot = cot / (4 * A); + T tmp(i, j, cot); + coef.push_back(tmp); + ii += cot; + } + } + T tmp(i, i, ii); + coef.push_back(tmp); + } + + L.setFromTriplets(coef.begin(), coef.end()); } diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp index 577e90f..ede0c2e 100644 --- a/src/massmatrix.cpp +++ b/src/massmatrix.cpp @@ -1,10 +1,48 @@ #include "massmatrix.h" -void massmatrix( - const Eigen::MatrixXd & l, - const Eigen::MatrixXi & F, - Eigen::DiagonalMatrix & M) -{ - // Add your code here +void massmatrix(const Eigen::MatrixXd & l, const Eigen::MatrixXi & F, + Eigen::DiagonalMatrix & M) { + // Add your code here + int vnum = F.maxCoeff() + 1; + int fnum = F.rows(); + + // calculate area of faces + Eigen::VectorXd FA(fnum); + for (int i = 0; i < fnum; i++) { + int v1 = F(i, 0); + int v2 = F(i, 1); + int v3 = F(i, 2); + double e1len = l(v1, v2); + double e2len = l(v1, v3); + double e3len = l(v2, v3); + double s = e1len + e2len + e3len; + s = s / 2; + double A2 = s * (s - e1len) * (s - e2len) * (s - e3len); + double A = sqrt(A2); + FA[i] = A; + } + + Eigen::MatrixXi V2F = Eigen::MatrixXi::Zero(vnum, fnum); + for (int i = 0; i < fnum; i++) { + int v1 = F(i, 0); + int v2 = F(i, 1); + int v3 = F(i, 2); + V2F(v1, i) = 1; + V2F(v2, i) = 1; + V2F(v3, i) = 1; + } + Eigen::VectorXd diag(vnum); + for (int i = 0; i < vnum; i++) { + double ii = 0; + for (int j = 0; j < fnum; j++) { + if (V2F(i, j) == 1) { + ii += FA[j]; + } + } + + diag[i] = ii / 3; + } + + M = diag.asDiagonal(); } diff --git a/src/smooth.cpp b/src/smooth.cpp index 6f72147..da74448 100644 --- a/src/smooth.cpp +++ b/src/smooth.cpp @@ -1,12 +1,35 @@ #include "smooth.h" +#include "massmatrix.h" +#include "cotmatrix.h" -void smooth( - const Eigen::MatrixXd & V, - const Eigen::MatrixXi & F, - const Eigen::MatrixXd & G, - double lambda, - Eigen::MatrixXd & U) -{ - // Replace with your code - U = G; +#include "igl/edge_lengths.h" +#include + +void smooth(const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, + const Eigen::MatrixXd & G, double lambda, Eigen::MatrixXd & U) { + // Replace with your code +// first, calculate face edges; + Eigen::MatrixXd l; + igl::edge_lengths(V, F, l); + + // second, calculate L + Eigen::SparseMatrix L; + cotmatrix(l, F, L); + + // third , M + Eigen::DiagonalMatrix M; + massmatrix(l, F, M); + + // solve equation + // next, we combine G with M + Eigen::SparseMatrix A = -lambda * L; + for (int i = 0; i < V.rows(); i++) { + A.coeffRef(i, i) += M.diagonal()(i); + } + // last, we solve Ax = B + Eigen::BiCGSTAB > solver; + solver.compute(A); + + Eigen::MatrixXd B = M * G; + U = solver.solve(B); }