From ccf589b32531e5999887bec950c0f38bf9e1d893 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Fri, 5 Apr 2024 23:27:14 -0400 Subject: add videos, concurrency & add dirty init --- README.md | 33 +++- readme-videos/armadillo.gif | Bin 1831687 -> 9152251 bytes readme-videos/bean.gif | Bin 2039976 -> 18613728 bytes readme-videos/bunny-conc.gif | Bin 0 -> 13272617 bytes readme-videos/bunny-norm.gif | Bin 0 -> 10736071 bytes readme-videos/peter.gif | Bin 6964791 -> 56931698 bytes readme-videos/sphere.gif | Bin 1258776 -> 7663441 bytes readme-videos/teapot.gif | Bin 1139394 -> 6681549 bytes readme-videos/tetrahedron.gif | Bin 478109 -> 5602789 bytes src/arap.cpp | 343 +++++++++++++++++++++++------------------- src/arap.h | 22 +++ 11 files changed, 244 insertions(+), 154 deletions(-) create mode 100644 readme-videos/bunny-conc.gif create mode 100644 readme-videos/bunny-norm.gif diff --git a/README.md b/README.md index 765c344..c637b9a 100644 --- a/README.md +++ b/README.md @@ -50,8 +50,27 @@ This sums to **95 points**. To score **100 points** (or more!), you’ll need to Your README **(5 points)** should describe: 1. How to run your program, such as how to load a specific mesh; and + +To run my program, you use the basic clion run configuration. Use release mode for best performance. + +There are two parameters you must modify: 1) the number of iterations and 2) the mesh to load. +Both are modifiable on the first two lines in the arap.cpp's init function. +Use an integer for the num iterations and a string to the path for the mesh. + 2. _Briefly_, all the features your code implements, including any known bugs. - - E.g.: "I implemented ARAP ... and affine progressive meshes ... however, my program doesn't work in these specific cases: ..." + +I implemented the basics of arap under the move() function. The order of the overall implementation is as follows: +- When the number of anchors changes: + - Estimate pprime from p. + - Compute the cotan weights. + - Define the constraints map. + - Compute the L matrix. +- For the number of iterations: + - Estimate the rotation matrices. + - Estimate the positions. +- Update the mesh vertices. + +The number of iterations is statically defined as a parameter. You should also have all the [example videos](#example-videos) described in the next section embedded in your README. @@ -67,3 +86,15 @@ For this project, we ask that you demonstrate to us that your program achieves t | You should be able to make the armadillo wave. | ![](./readme-videos/armadillo.gif) | | Attempting to deform `tetrahedron.obj` should not cause it to collapse or behave erratically. | ![](./readme-videos/tetrahedron.gif) | | Attempting to deform a (large) mesh like `bunny.obj` or `peter.obj` should not cause your code to crash. | ![](./readme-videos/peter.gif) | + +### Extra Feature: Parallelization + +For the solving algorithm, I divided the dimensions and parallelized my code using the QTFramework to concurrently calculate the estimated positions. +See line 280 in arap.c for the implementation. +The video below shows both implementations (concurrent and normal) when deforming bunny.obj, a large mesh, at 1000 iterations. +It appears the concurrent implementation is faster by only about a second. + +| Program Specification | Video | +|:------------------------------------------------------------------|:------------------------------------| +| Modifying Peter at 1000 iterations with no concurrent framework. | ![](./readme-videos/bunny-conc.gif) | +| Modifying Peter at 1000 iterations using QTConcurrent map-reduce. | ![](./readme-videos/bunny-norm.gif) | diff --git a/readme-videos/armadillo.gif b/readme-videos/armadillo.gif index 5b15b34..26c75ed 100644 Binary files a/readme-videos/armadillo.gif and b/readme-videos/armadillo.gif differ diff --git a/readme-videos/bean.gif b/readme-videos/bean.gif index 5a80c0b..be92079 100644 Binary files a/readme-videos/bean.gif and b/readme-videos/bean.gif differ diff --git a/readme-videos/bunny-conc.gif b/readme-videos/bunny-conc.gif new file mode 100644 index 0000000..459c927 Binary files /dev/null and b/readme-videos/bunny-conc.gif differ diff --git a/readme-videos/bunny-norm.gif b/readme-videos/bunny-norm.gif new file mode 100644 index 0000000..4774bbb Binary files /dev/null and b/readme-videos/bunny-norm.gif differ diff --git a/readme-videos/peter.gif b/readme-videos/peter.gif index bb00f4b..daa6ce2 100644 Binary files a/readme-videos/peter.gif and b/readme-videos/peter.gif differ diff --git a/readme-videos/sphere.gif b/readme-videos/sphere.gif index 47a131a..e0a8ae1 100644 Binary files a/readme-videos/sphere.gif and b/readme-videos/sphere.gif differ diff --git a/readme-videos/teapot.gif b/readme-videos/teapot.gif index 7c202bb..e89fcba 100644 Binary files a/readme-videos/teapot.gif and b/readme-videos/teapot.gif differ diff --git a/readme-videos/tetrahedron.gif b/readme-videos/tetrahedron.gif index 866f1e2..0a8f139 100644 Binary files a/readme-videos/tetrahedron.gif and b/readme-videos/tetrahedron.gif differ diff --git a/src/arap.cpp b/src/arap.cpp index 4ba3bc8..f4ac929 100644 --- a/src/arap.cpp +++ b/src/arap.cpp @@ -6,10 +6,8 @@ #include #include -#include -#include -#include -#include +#include +#include using namespace std; using namespace Eigen; @@ -18,11 +16,15 @@ ARAP::ARAP() {} void ARAP::init(Eigen::Vector3f &coeffMin, Eigen::Vector3f &coeffMax) { + // PARAMETERS SET HERE + m_num_iterations = 5; + m_mesh_path = "meshes/bean.obj"; + vector vertices; vector triangles; // If this doesn't work for you, remember to change your working directory - if (MeshLoader::loadTriMesh("meshes/sphere.obj", vertices, triangles)) { + if (MeshLoader::loadTriMesh(m_mesh_path, vertices, triangles)) { m_shape.init(vertices, triangles); } @@ -34,9 +36,17 @@ void ARAP::init(Eigen::Vector3f &coeffMin, Eigen::Vector3f &coeffMax) } coeffMin = all_vertices.colwise().minCoeff(); coeffMax = all_vertices.colwise().maxCoeff(); -} + // initialize the number of anchors + num_anchors = m_shape.getAnchors().size(); + // initialze DS for ARAP + m_rotations = vector(m_shape.getVertices().size(), RM::Identity(3, 3)); + m_vtx_to_free_vtx = vector(m_shape.getVertices().size(), -1); + m_edge_weights = Sparse(m_shape.getVertices().size(), m_shape.getVertices().size()); + m_b = PM::Zero(3, m_shape.getVertices().size()); + m_b_fixed = PM::Zero(3, m_shape.getVertices().size()); +} // Move an anchored vertex, defined by its index, to targetPosition void ARAP::move(int vertex, Vector3f targetPosition) @@ -44,14 +54,11 @@ void ARAP::move(int vertex, Vector3f targetPosition) std::vector new_vertices = m_shape.getVertices(); const std::unordered_set& anchors = m_shape.getAnchors(); - std::cout << "anchors size" << anchors.size() << std::endl; - std::cout << "targetPosition" << targetPosition << std::endl; - - // TODO: implement ARAP here - new_vertices[vertex] = targetPosition; + int new_num_anchors = anchors.size(); + // TODO: implement ARAP here + new_vertices[vertex] = targetPosition; // make a copy of the vertices into the matrix points - typedef Matrix PM; // position matrix PM p(3, m_shape.getVertices().size()); for (int i = 0; i < m_shape.getVertices().size(); ++i) { @@ -64,153 +71,158 @@ void ARAP::move(int vertex, Vector3f targetPosition) p_prime.col(i) = new_vertices[i]; } - // compute the cotan weights - typedef Triplet Tri; - vector triplets; - triplets.reserve(3 * m_shape.getFaces().size()); - for (DenseIndex f = 0; f < m_shape.getFaces().size(); ++f) - { - Vector3i face = m_shape.getFaces()[f]; - Vector3f v0 = m_shape.getVertices()[face[0]]; - Vector3f v1 = m_shape.getVertices()[face[1]]; - Vector3f v2 = m_shape.getVertices()[face[2]]; - - Vector3f e0 = v1 - v0; - Vector3f e1 = v2 - v1; - Vector3f e2 = v0 - v2; - - // protect against degen triangles - float s0 = e0.squaredNorm(); - s0 = max(s0, 1e-8f); - float s1 = e1.squaredNorm(); - s1 = max(s1, 1e-8f); - float s2 = e2.squaredNorm(); - s2 = max(s2, 1e-8f); - - s0 = sqrt(s0); - s1 = sqrt(s1); - s2 = sqrt(s2); - - float semi_perimeter = (s0 + s1 + s2) * 0.5f; - float area = sqrt(semi_perimeter * (semi_perimeter - s0) * (semi_perimeter - s1) * (semi_perimeter - s2)); - area = max(area, 1e-8f); - - float cot0 = (-s0 * s0 + s1 * s1 + s2 * s2) / (4.0f * area); - cot0 = max(cot0, 1e-10f) * 0.5f; - float cot1 = (s0 * s0 - s1 * s1 + s2 * s2) / (4.0f * area); - cot1 = max(cot1, 1e-10f) * 0.5f; - float cot2 = (s0 * s0 + s1 * s1 - s2 * s2) / (4.0f * area); - cot2 = max(cot2, 1e-10f) * 0.5f; - - pair edge0 = - face[0] > face[1] ? make_pair(face[1], face[0]) : make_pair(face[0], face[1]); - pair edge1 = - face[1] > face[2] ? make_pair(face[2], face[1]) : make_pair(face[1], face[2]); - pair edge2 = - face[2] > face[0] ? make_pair(face[0], face[2]) : make_pair(face[2], face[0]); - - triplets.push_back(Tri(edge0.first, edge0.second, cot0)); - triplets.push_back(Tri(edge1.first, edge1.second, cot1)); - triplets.push_back(Tri(edge2.first, edge2.second, cot2)); - triplets.push_back(Tri(edge0.second, edge0.first, cot0)); - triplets.push_back(Tri(edge1.second, edge1.first, cot1)); - triplets.push_back(Tri(edge2.second, edge2.first, cot2)); - } - // build the sparse matrix of edge weights - SparseMatrix edge_weights(m_shape.getVertices().size(), m_shape.getVertices().size()); - edge_weights.reserve(VectorXi::Constant(m_shape.getVertices().size(), 7)); - edge_weights.setZero(); - edge_weights.setFromTriplets(triplets.begin(), triplets.end()); - - // initialize the rotation matrices - typedef Matrix R; - // TODO: move this to only init once - vector rotations(m_shape.getVertices().size(), R::Identity(3, 3)); -// rotations.clear(); -// rotations.resize(m_shape.getVertices().size(), R::Identity(3, 3)); - - // make a map that reduces the constrained vertices (vtx -> free vtx) - vector vtx_to_free_vtx = vector(m_shape.getVertices().size(), -1); - // vtx_to_free_vtx.reserve(m_shape.getVertices().size()); - Index count_free = 0; - for (int i = 0; i < m_shape.getVertices().size(); ++i) + if (new_num_anchors != num_anchors) { - if (!anchors.contains(i)) + // compute the cotan weights + typedef Triplet Tri; + vector triplets; + triplets.reserve(3 * m_shape.getFaces().size()); + for (DenseIndex f = 0; f < m_shape.getFaces().size(); ++f) { - // if were not an anchor, we are free - vtx_to_free_vtx[i] = count_free; - count_free++; - } - } - - // update pprime with the constraints into the linear system - for (Index i = 0; i < m_shape.getVertices().size(); ++i) - { - if (vtx_to_free_vtx[i] == -1) - { - // if we're an anchor, set this position - p_prime.col(i) = new_vertices[i]; - } - } + Vector3i face = m_shape.getFaces()[f]; + Vector3f v0 = m_shape.getVertices()[face[0]]; + Vector3f v1 = m_shape.getVertices()[face[1]]; + Vector3f v2 = m_shape.getVertices()[face[2]]; + + Vector3f e0 = v1 - v0; + Vector3f e1 = v2 - v1; + Vector3f e2 = v0 - v2; + + // protect against degen triangles + float s0 = e0.squaredNorm(); + s0 = max(s0, 1e-8f); + float s1 = e1.squaredNorm(); + s1 = max(s1, 1e-8f); + float s2 = e2.squaredNorm(); + s2 = max(s2, 1e-8f); + + s0 = sqrt(s0); + s1 = sqrt(s1); + s2 = sqrt(s2); + + float semi_perimeter = (s0 + s1 + s2) * 0.5f; + float area = sqrt(semi_perimeter * (semi_perimeter - s0) * (semi_perimeter - s1) * (semi_perimeter - s2)); + area = max(area, 1e-8f); + + float cot0 = (-s0 * s0 + s1 * s1 + s2 * s2) / (4.0f * area); + cot0 = max(cot0, 1e-10f) * 0.5f; + float cot1 = (s0 * s0 - s1 * s1 + s2 * s2) / (4.0f * area); + cot1 = max(cot1, 1e-10f) * 0.5f; + float cot2 = (s0 * s0 + s1 * s1 - s2 * s2) / (4.0f * area); + cot2 = max(cot2, 1e-10f) * 0.5f; + + pair edge0 = + face[0] > face[1] ? make_pair(face[1], face[0]) : make_pair(face[0], face[1]); + pair edge1 = + face[1] > face[2] ? make_pair(face[2], face[1]) : make_pair(face[1], face[2]); + pair edge2 = + face[2] > face[0] ? make_pair(face[0], face[2]) : make_pair(face[2], face[0]); + + triplets.push_back(Tri(edge0.first, edge0.second, cot0)); + triplets.push_back(Tri(edge1.first, edge1.second, cot1)); + triplets.push_back(Tri(edge2.first, edge2.second, cot2)); + triplets.push_back(Tri(edge0.second, edge0.first, cot0)); + triplets.push_back(Tri(edge1.second, edge1.first, cot1)); + triplets.push_back(Tri(edge2.second, edge2.first, cot2)); + + // build the sparse matrix of edge weights + m_edge_weights.resize(m_shape.getVertices().size(), m_shape.getVertices().size()); + m_edge_weights.reserve(VectorXi::Constant(m_shape.getVertices().size(), 7)); + m_edge_weights.setZero(); + m_edge_weights.setFromTriplets(triplets.begin(), triplets.end()); + + // initialize the rotation matrices + // TODO: move this to only init once +// vector rotations(m_shape.getVertices().size(), RM::Identity(3, 3)); +// rotations.clear(); +// rotations.resize(m_shape.getVertices().size(), RM::Identity(3, 3)); + + // make a map that reduces the constrained vertices (vtx -> free vtx) + m_vtx_to_free_vtx.clear(); + m_vtx_to_free_vtx.reserve(m_shape.getVertices().size()); + Index count_free = 0; + for (int i = 0; i < m_shape.getVertices().size(); ++i) + { + if (!anchors.contains(i)) + { + // if were not an anchor, we are free + m_vtx_to_free_vtx[i] = count_free; + count_free++; + } + } - // setup the linear system we will use to solve - SparseMatrix L(count_free, count_free); - // L.resize(count_free, count_free); - L.reserve(VectorXi::Constant(count_free, 7)); - L.setZero(); + // update pprime with the constraints into the linear system + for (Index i = 0; i < m_shape.getVertices().size(); ++i) + { + if (m_vtx_to_free_vtx[i] == -1) + { + // if we're an anchor, set this position + p_prime.col(i) = new_vertices[i]; + } + } - // setup bfixed - PM b_fixed = PM::Zero(3, count_free); - // b_fixed.resize(3, count_free); + // setup the linear system we will use to solve + SparseMatrix L(count_free, count_free); + // L.resize(count_free, count_free); + L.reserve(VectorXi::Constant(count_free, 7)); + L.setZero(); - // make the L matrix using triplet method - vector triplets_L; - triplets_L.reserve(7 * count_free); - for (Index i = 0; i < edge_weights.outerSize(); i++) - { - Index free_index = vtx_to_free_vtx[i]; - if (free_index == -1) - { - // do not add constraints to the linear system - continue; - } + // setup bfixed + m_b_fixed = PM::Zero(3, count_free); + m_b_fixed.resize(3, count_free); - for (SparseMatrix::InnerIterator it(edge_weights, i); it; ++it) - { - Index j = it.col(); - Index free_index_j = vtx_to_free_vtx[j]; - float wij = it.value(); - if (free_index_j == -1) + // make the L matrix using triplet method + vector triplets_L; + triplets_L.reserve(7 * count_free); + for (Index i = 0; i < m_edge_weights.outerSize(); i++) { - // if the neighbor is an anchor, add to b_fixed - Vector3f location = new_vertices[j]; - b_fixed.col(free_index) += wij * location; - } - else - { - triplets_L.push_back(Tri(free_index, free_index_j, -wij)); + Index free_index = m_vtx_to_free_vtx[i]; + if (free_index == -1) + { + // do not add constraints to the linear system + continue; + } + + for (SparseMatrix::InnerIterator it(m_edge_weights, i); it; ++it) + { + Index j = it.col(); + Index free_index_j = m_vtx_to_free_vtx[j]; + float wij = it.value(); + if (free_index_j == -1) + { + // if the neighbor is an anchor, add to b_fixed + Vector3f location = new_vertices[j]; + m_b_fixed.col(free_index) += wij * location; + } + else + { + triplets_L.push_back(Tri(free_index, free_index_j, -wij)); + } + triplets_L.push_back(Tri(free_index, free_index, wij)); + } } - triplets_L.push_back(Tri(free_index, free_index, wij)); + // create and solve the L matrix + L.setFromTriplets(triplets_L.begin(), triplets_L.end()); + m_solver.compute(L); } + } - // create and solve the L matrix - L.setFromTriplets(triplets_L.begin(), triplets_L.end()); - SimplicialLDLT> solver; - solver.compute(L); - for (int iter = 0; iter < 5; iter++) + // do the iterations + for (int iter = 0; iter < m_num_iterations; iter++) { // estimate the rotations - rotations.clear(); - rotations.resize(m_shape.getVertices().size(), R::Identity(3, 3)); + m_rotations.clear(); + m_rotations.resize(m_shape.getVertices().size(), RM::Identity(3, 3)); + typedef Matrix S; - for (Index i = 0; i < edge_weights.outerSize(); ++i) + for (Index i = 0; i < m_edge_weights.outerSize(); ++i) { const auto &pi = p.col(i); const auto &ppi = p_prime.col(i); S cov = S::Zero(3, 3); - for (SparseMatrix::InnerIterator it(edge_weights, i); it; ++it) + for (SparseMatrix::InnerIterator it(m_edge_weights, i); it; ++it) { // get the weight float wij = it.value(); @@ -231,50 +243,75 @@ void ARAP::move(int vertex, Vector3f targetPosition) S I = S::Identity(3, 3); I(2, 2) = (v * u_trans).determinant(); - rotations[i] = v * I * u_trans; + m_rotations[i] = v * I * u_trans; } // estimate the positions // make a copy of b_fixed to modify - PM b = b_fixed; - for (Index i = 0; i < edge_weights.outerSize(); ++i) + m_b = m_b_fixed; + for (Index i = 0; i < m_edge_weights.outerSize(); ++i) { // TODO add constraint check here - Index free_index = vtx_to_free_vtx[i]; + Index free_index = m_vtx_to_free_vtx[i]; if (free_index == -1) { // do not consider constraints continue; } - for (SparseMatrix::InnerIterator it(edge_weights, i); it; ++it) + for (SparseMatrix::InnerIterator it(m_edge_weights, i); it; ++it) { // get the weight float wij = it.value(); // get the rotation matrices for the vertices - const S &R = rotations[i]; - const S &Rj = rotations[it.col()]; + const S &R = m_rotations[i]; + const S &Rj = m_rotations[it.col()]; // calculate the energy Eigen::Matrix pt = (R + Rj) * (p.col(i) - p.col(it.col())); pt *= wij * 0.5f; - b.col(free_index) += pt; + m_b.col(free_index) += pt; } } // solve the linear system for each dimension - Matrix U; +// Matrix U; +// for (int j = 0; j < 3; ++j) +// { +// U = solver.solve(b.row(j).transpose()); +// +// Index idx = 0; +// for (int k = 0; k < m_shape.getVertices().size(); ++k) +// { +// if (vtx_to_free_vtx[k] != -1) +// { +// p_prime(j, k) = U(idx); +// idx++; +// } +// } +// } + +// // EXTRA FEATURE - parallelize the solving of the linear system by dimension + QList dims {0, 1, 2}; + QList> results = QtConcurrent::blockingMapped( + dims, + [&](int j) + { + Matrix U; + U = m_solver.solve(m_b.row(j).transpose()); + return U; + }); + + // fill the results into p_prime for (int j = 0; j < 3; ++j) { - U = solver.solve(b.row(j).transpose()); - Index idx = 0; for (int k = 0; k < m_shape.getVertices().size(); ++k) { - if (vtx_to_free_vtx[k] != -1) + if (m_vtx_to_free_vtx[k] != -1) { - p_prime(j, k) = U(idx); + p_prime(j, k) = results[j](idx); idx++; } } diff --git a/src/arap.h b/src/arap.h index afa7d63..10a2904 100644 --- a/src/arap.h +++ b/src/arap.h @@ -3,6 +3,12 @@ #include "graphics/shape.h" #include "Eigen/StdList" #include "Eigen/StdVector" +#include +#include +#include +#include +#include +#include class Shader; @@ -43,4 +49,20 @@ public: { return m_shape.getAnchorPos(lastSelected, pos, ray, start); } + + // for determinig when to recompute + int num_anchors; + + typedef Eigen::Matrix PM; // position matrix + typedef Eigen::Matrix RM; // rotation matrix + typedef Eigen::SparseMatrix Sparse; + + std::vector m_rotations; + Sparse m_edge_weights; + PM m_b, m_b_fixed; + std::vector m_vtx_to_free_vtx; + Eigen::SimplicialLDLT> m_solver; + + int m_num_iterations; + const char * m_mesh_path; }; -- cgit v1.2.3-70-g09d2