From c7b575b513235bea488d3d9720ecab11dcd21c3b Mon Sep 17 00:00:00 2001 From: andrew Date: Tue, 13 Feb 2018 01:22:28 -0500 Subject: [PATCH] Finished the assignment --- src/cotmatrix.cpp | 75 ++++++++++++++++++++++++++++++++++++++++++---- src/massmatrix.cpp | 35 ++++++++++++++++++---- src/smooth.cpp | 43 +++++++++++++++++++++----- 3 files changed, 135 insertions(+), 18 deletions(-) diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp index c888e5c..08daf24 100644 --- a/src/cotmatrix.cpp +++ b/src/cotmatrix.cpp @@ -1,10 +1,75 @@ #include "cotmatrix.h" +/** +* @brief Solves with the cosine law to produce the cosine of the angle +* opposite to the edge c. +*/ +double cosineLaw(double a, double b, double c) { + double top = a * a + b * b - c * c; + double bottom = 2 * a * b; + return top / bottom; +} + +/** +* @brief Computes the cotagenent of the angle opposite to the edges with the +* given lengths in a triangle. +* +* @param lengths input array of 3 edge lengths. +* @param angles output array of 3 angles. +*/ +void lengthsToCotangents(double* lengths, double* cots) { + + // First we apply cosine law. + double cosines[3]; + cosines[0] = cosineLaw(lengths[1], lengths[2], lengths[0]); + cosines[1] = cosineLaw(lengths[2], lengths[0], lengths[1]); + cosines[2] = cosineLaw(lengths[0], lengths[1], lengths[2]); + + // Now we apply pythagorean and cotangent identities to produce the + // cotangent. + for (int i = 0; i < 3; i++) { + + // Pythaogrean identity is valid because the angle is known to be + // less than PI; so sine wouldn't have gone negative in there. + double cosine = cosines[i]; + double sine = sqrt(1 - cosine * cosine); + double cotangent = cosine / sine; + cots[i] = cotangent; + } +} + void cotmatrix( - const Eigen::MatrixXd & l, - const Eigen::MatrixXi & F, - Eigen::SparseMatrix & L) -{ - // Add your code here + const Eigen::MatrixXd & l, + const Eigen::MatrixXi & F, + Eigen::SparseMatrix & L) { + + int vertexCount = F.maxCoeff() + 1; + int faceCount = l.rows(); + + L.resize(vertexCount, vertexCount); + for (int i = 0; i < faceCount; i++) { + + // Find cotangents opposite to each edge. + double lengths[] = { l(i, 0), l(i, 1), l(i, 2) }; + double cots[3]; + lengthsToCotangents(lengths, cots); + + // For each edge, fill in the appropriate matrix stuff. + L.coeffRef(F(i, 0), F(i, 1)) += 0.5 * cots[2]; + L.coeffRef(F(i, 1), F(i, 2)) += 0.5 * cots[0]; + L.coeffRef(F(i, 2), F(i, 0)) += 0.5 * cots[1]; + + // And its symmetric after all... + L.coeffRef(F(i, 1), F(i, 0)) += 0.5 * cots[2]; + L.coeffRef(F(i, 2), F(i, 1)) += 0.5 * cots[0]; + L.coeffRef(F(i, 0), F(i, 2)) += 0.5 * cots[1]; + } + + // Now set the diagonal + for (int i = 0; i < vertexCount; i++) { + L.insert(i, i) = -L.row(i).sum(); + } + + } diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp index 577e90f..bb6b887 100644 --- a/src/massmatrix.cpp +++ b/src/massmatrix.cpp @@ -1,10 +1,35 @@ #include "massmatrix.h" void massmatrix( - const Eigen::MatrixXd & l, - const Eigen::MatrixXi & F, - Eigen::DiagonalMatrix & M) -{ - // Add your code here + const Eigen::MatrixXd & l, + const Eigen::MatrixXi & F, + Eigen::DiagonalMatrix & M) { + + int vertexCount = F.maxCoeff() + 1; + int triangleCount = F.rows(); + + M.resize(vertexCount); + for (int i = 0; i < vertexCount; i++) { + + double totalArea = 0; + for (int t = 0; t < triangleCount; t++) { + + // Check if the vertex i is contained in the triangle t. + if (F(t, 0) != i && F(t, 1) != i && F(t, 2) != i) { + continue; + } + + // Compute area of the triangle from the side lengths. + // Use Heron's formula. + double s = (l(i, 0) + l(i, 1) + l(i, 2)) / 2; + double area = sqrt(s * (s - l(i, 0)) * (s - l(i, 1)) * + (s - l(i, 2))); + totalArea += area; + } + M.diagonal()[i] = totalArea; + } + + M = M * (1 / 3.0); + } diff --git a/src/smooth.cpp b/src/smooth.cpp index 6f72147..77ca73c 100644 --- a/src/smooth.cpp +++ b/src/smooth.cpp @@ -1,12 +1,39 @@ #include "smooth.h" +#include "igl/edge_lengths.h" +#include "cotmatrix.h" +#include "massmatrix.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; + const Eigen::MatrixXd & V, + const Eigen::MatrixXi & F, + const Eigen::MatrixXd & G, + double lambda, + Eigen::MatrixXd & U) { + + // Compute matrix of edge lengths. + Eigen::MatrixXd edgeLengths; + edgeLengths.resizeLike(F); + igl::edge_lengths(V, F, edgeLengths); + + // Compute Laplacian and Mass Matrices + Eigen::SparseMatrix laplacian; + Eigen::DiagonalMatrix mass; + cotmatrix(edgeLengths, F, laplacian); + massmatrix(edgeLengths, F, mass); + + // Treat the system given as a linear system we can solve with Ax = b. + // I construct A weirdly because of the data-types involved. + Eigen::MatrixXd b = mass * G; + Eigen::SparseMatrix A = -lambda * laplacian; + for (int i = 0; i < mass.rows(); i++) { + A.coeffRef(i, i) += mass.diagonal()[i]; + } + + // Now solve A to obtain G for the next iteration. + Eigen::ConjugateGradient> solver; + solver.compute(A); + + U.resizeLike(G); + U = solver.solve(b); + return; }