From c2b474c81112f665fb552d3aa9d2e7ab7b5b02b4 Mon Sep 17 00:00:00 2001 From: Mark Date: Tue, 13 Nov 2018 01:27:43 -0500 Subject: [PATCH 1/3] mean curvature --- src/mean_curvature.cpp | 33 +++++++++++++++++++++++++++++++-- src/principal_curvatures.cpp | 10 +++++----- 2 files changed, 36 insertions(+), 7 deletions(-) diff --git a/src/mean_curvature.cpp b/src/mean_curvature.cpp index 41fce89..8b35e24 100644 --- a/src/mean_curvature.cpp +++ b/src/mean_curvature.cpp @@ -1,10 +1,39 @@ #include "../include/mean_curvature.h" +#include +#include +#include + void mean_curvature( const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & H) { - // Replace with your code - H = Eigen::VectorXd::Zero(V.rows()); + // Compute mean curvature normal `Hn = M^{-1}*L*V` + // by solving M*Hn = L*V + + Eigen::SparseMatrix L, M; + igl::cotmatrix(V, F, L); + igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M); + + Eigen::MatrixXd b; + b = L*V; + Eigen::MatrixXd Hn(V.rows(), 3); + Eigen::SimplicialLDLT> solver; + solver.compute(M); + Hn = solver.solve(b); + + // Determine sign + H.resize(V.rows()); + + Eigen::MatrixXd N; + igl::per_vertex_normals(V, F, N); + + for (int i = 0; i < V.rows(); ++i) { + H(i) = Hn.row(i).norm(); + if (Hn.row(i).dot(N.row(i)) > 0) { + H(i) = -H(i); + } + } } + diff --git a/src/principal_curvatures.cpp b/src/principal_curvatures.cpp index 6a51d2e..a750a56 100644 --- a/src/principal_curvatures.cpp +++ b/src/principal_curvatures.cpp @@ -8,9 +8,9 @@ void principal_curvatures( Eigen::VectorXd & K1, Eigen::VectorXd & K2) { - // Replace with your code - K1 = Eigen::VectorXd::Zero(V.rows()); - K2 = Eigen::VectorXd::Zero(V.rows()); - D1 = Eigen::MatrixXd::Zero(V.rows(),3); - D2 = Eigen::MatrixXd::Zero(V.rows(),3); + // Replace with your code + K1 = Eigen::VectorXd::Zero(V.rows()); + K2 = Eigen::VectorXd::Zero(V.rows()); + D1 = Eigen::MatrixXd::Zero(V.rows(),3); + D2 = Eigen::MatrixXd::Zero(V.rows(),3); } From 77f3b1bfe95154ea62d47e5d81244b1f0dd5f50e Mon Sep 17 00:00:00 2001 From: Mark Date: Tue, 13 Nov 2018 02:02:42 -0500 Subject: [PATCH 2/3] angle defect --- src/angle_defect.cpp | 24 ++++++++++++++++++++++-- src/internal_angles.cpp | 25 +++++++++++++++++++++++-- 2 files changed, 45 insertions(+), 4 deletions(-) diff --git a/src/angle_defect.cpp b/src/angle_defect.cpp index 78589fb..3e4519b 100644 --- a/src/angle_defect.cpp +++ b/src/angle_defect.cpp @@ -1,9 +1,29 @@ #include "../include/angle_defect.h" +#include "../include/internal_angles.h" + +#include void angle_defect( const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & D) { - D = Eigen::VectorXd::Zero(V.rows()); -} + // Interior angles + + Eigen::MatrixXd l_sqr, A; + igl::squared_edge_lengths(V, F, l_sqr); + internal_angles(l_sqr, A); + + // Compute angle defect + // \Sigma_{i} = 2\pi - \sum_{f\in F(i)} InteriorAngle(i, f) + + D = Eigen::VectorXd::Ones(V.rows()); + D = 2 * M_PI * D; + + for (int i = 0; i < F.rows(); ++i) { + for (int j = 0; j < 3; ++j) { + auto v = F(i, j); + D(v) = D(v) - A(i, j); + } + } +} \ No newline at end of file diff --git a/src/internal_angles.cpp b/src/internal_angles.cpp index 3e14c75..511c75e 100644 --- a/src/internal_angles.cpp +++ b/src/internal_angles.cpp @@ -1,8 +1,29 @@ #include "../include/internal_angles.h" +#include + void internal_angles( const Eigen::MatrixXd & l_sqr, Eigen::MatrixXd & A) { - // Add with your code -} + A.resizeLike(l_sqr); + + // return angle opposite of edge `c` + auto angle = [](double a2, double b2, double c2) { + return std::acos((a2 + b2 - c2) / (2 * sqrt(a2) * sqrt(b2))); + }; + + double a2, b2, c2; + for (int i = 0; i < l_sqr.rows(); ++i) + { + // row = [A, B, C] where `A` opposite of edge `a` + + a2 = l_sqr(i, 0); + b2 = l_sqr(i, 1); + c2 = l_sqr(i, 2); + + A(i, 0) = angle(b2, c2, a2); + A(i, 1) = angle(c2, a2, b2); + A(i, 2) = angle(a2, b2, c2); + } +} \ No newline at end of file From 6aac4f9e93786d876329ff20718733d5afcaea72 Mon Sep 17 00:00:00 2001 From: Mark Date: Tue, 13 Nov 2018 19:55:21 -0500 Subject: [PATCH 3/3] add --- .vscode/settings.json | 19 ++++++ CMakeLists.txt | 4 ++ src/principal_curvatures.cpp | 108 ++++++++++++++++++++++++++++++++++- 3 files changed, 130 insertions(+), 1 deletion(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..c637a9f --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,19 @@ +{ + "files.associations": { + "*.txt": "cmake", + "complex": "cpp", + "__config": "cpp", + "__nullptr": "cpp", + "cstddef": "cpp", + "exception": "cpp", + "initializer_list": "cpp", + "new": "cpp", + "stdexcept": "cpp", + "type_traits": "cpp", + "typeinfo": "cpp", + "algorithm": "cpp", + "hashtable": "cpp", + "unordered_map": "cpp", + "utility": "cpp" + } +} \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 76d0c5b..0ec3e57 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,8 @@ cmake_minimum_required(VERSION 2.6) project(curvature) + +add_definitions(${CMAKE_CXX_FLAGS} "-g") +add_definitions(${CMAKE_CXX_FLAGS} "-Wall") + set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/shared/cmake) include("${CMAKE_CURRENT_SOURCE_DIR}/shared/cmake/CMakeLists.txt") diff --git a/src/principal_curvatures.cpp b/src/principal_curvatures.cpp index a750a56..6bacc46 100644 --- a/src/principal_curvatures.cpp +++ b/src/principal_curvatures.cpp @@ -1,4 +1,10 @@ #include "../include/principal_curvatures.h" +#include + +#include +#include + +#include void principal_curvatures( const Eigen::MatrixXd & V, @@ -8,9 +14,109 @@ void principal_curvatures( Eigen::VectorXd & K1, Eigen::VectorXd & K2) { - // Replace with your code K1 = Eigen::VectorXd::Zero(V.rows()); K2 = Eigen::VectorXd::Zero(V.rows()); D1 = Eigen::MatrixXd::Zero(V.rows(),3); D2 = Eigen::MatrixXd::Zero(V.rows(),3); + + + // Gather two-ring neighborhood of each vertex + + Eigen::SparseMatrix A; + igl::adjacency_matrix(F, A); + + // Compute 2-ring by matrix multiplication + // vertices in 1-ring of v also in 2-ring of v by nature of triangle mesh + // B(i,i) > 0, since every vertex is distance 2 from itself. + Eigen::SparseMatrix B; + B = A*A; + + // Construct neighborhood distance matrix + int k = 0; + Eigen::MatrixXd P; + + // Eigen Decomposition solver + Eigen::SelfAdjointEigenSolver solver1; + // Eigen-vectors and projected positions + Eigen::MatrixXd W, PW; + + // Fitting height surface with quadratic polynomial + Eigen::MatrixXd U, Uinv; + Eigen::VectorXd x(5); + + // Shape operator `S`, foundamental forms `ff1`, `ff2` + Eigen::Matrix2d ff1, ff2, S; + double e,f,g,E,FF,G; + Eigen::SelfAdjointEigenSolver solver2; + + // principle tangent + Eigen::Vector2d pt1, pt2; + + + for (int i = 0; i < V.rows(); ++i) { + + // Construct `P`, distance from `v` for vertices in 2-ring of `v` + + P.resize(B.innerVector(i).nonZeros(), 3); + k = 0; + for (Eigen::SparseMatrix::InnerIterator it(B, i); it; ++it) { + // Visit nonzero element in i-th column + P.row(k) = V.row(it.row()) - V.row(i); + k += 1; + } + + // PCA on `P` + + solver1.compute(P.transpose() * P); + W = solver1.eigenvectors(); + PW = P * W; + + // Note `W` corresponds to eigenvalue in increasing order + // hence first principle component is W.col(2), i.e. u + + // Fit polynomial to height-field surface + // by solving a least squared problem for w=a₁u+a₂v+a₃u²+a₄uv+a₅v² + // i.e. find x = [a₁ a₂ a₃ a₄ a₅]^T by solving `Ux = b` where + // `U` consists of value of `u,v,u²,uv,v²` and `b = PW.col(0)` + + U.resize(P.rows(), 5); + U.col(0) = PW.col(2); + U.col(1) = PW.col(1); + U.col(2) = PW.col(2).array().square(); + U.col(3) = PW.col(2).array() * PW.col(1).array(); + U.col(4) = PW.col(1).array().square(); + + igl::pinv(U, Uinv); + x = Uinv * PW.col(0); + + // Compute shape operator `S` + + E = 1 + x(0)*x(0); + FF = x(0)*x(1); + G = 1 + x(1)*x(1); + e = 2*x(2) / sqrt(E + G - 1); + f = x(3) / sqrt(E + G - 1); + g = 2*x(4) / sqrt(E + G - 1); + + ff2 << e, f, + f, g; + ff1 << E, FF, + FF, G; + S = -ff2 * ff1.inverse(); + + // Eigen-decomposition on S + // eigenvalues are principle curvatures + // eigenvectors are principle tangent directions in (u,v) basis + solver2.compute(S.transpose() * S); + + K1(i) = solver2.eigenvalues()(1); + K2(i) = solver2.eigenvalues()(0); + + // map pt = (a,b) \in span{u,v} to 3D: au + bv + pt1 = solver2.eigenvectors().col(1); + pt2 = solver2.eigenvectors().col(0); + + D1.row(i) = pt1(0) * W.col(2) + pt1(1) * W.col(1); + D2.row(i) = pt2(0) * W.col(2) + pt2(1) * W.col(1); + } }