aboutsummaryrefslogtreecommitdiff
path: root/src/arap.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/arap.cpp')
-rw-r--r--src/arap.cpp343
1 files changed, 190 insertions, 153 deletions
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 <map>
#include <vector>
-#include <Eigen/Core>
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-#include <Eigen/SVD>
+#include <QList>
+#include <QtConcurrent>
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<Vector3f> vertices;
vector<Vector3i> 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<RM>(m_shape.getVertices().size(), RM::Identity(3, 3));
+ m_vtx_to_free_vtx = vector<Index>(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<Eigen::Vector3f> new_vertices = m_shape.getVertices();
const std::unordered_set<int>& 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<float, 3, Dynamic> 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<float> Tri;
- vector<Tri> 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<int, int> edge0 =
- face[0] > face[1] ? make_pair(face[1], face[0]) : make_pair(face[0], face[1]);
- pair<int, int> edge1 =
- face[1] > face[2] ? make_pair(face[2], face[1]) : make_pair(face[1], face[2]);
- pair<int, int> 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<float, RowMajor> 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<float, 3, 3> R;
- // TODO: move this to only init once
- vector<R> 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<Index> vtx_to_free_vtx = vector<Index>(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<float> Tri;
+ vector<Tri> 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<int, int> edge0 =
+ face[0] > face[1] ? make_pair(face[1], face[0]) : make_pair(face[0], face[1]);
+ pair<int, int> edge1 =
+ face[1] > face[2] ? make_pair(face[2], face[1]) : make_pair(face[1], face[2]);
+ pair<int, int> 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<RM> 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<float, RowMajor> 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<float, RowMajor> 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<Tri> 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<float, RowMajor>::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<Tri> 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<float, RowMajor>::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<SparseMatrix<float>> 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<float, 3, 3> 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<float, RowMajor>::InnerIterator it(edge_weights, i); it; ++it)
+ for (SparseMatrix<float, RowMajor>::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<float, RowMajor>::InnerIterator it(edge_weights, i); it; ++it)
+ for (SparseMatrix<float, RowMajor>::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<float, 3, 1> 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<float, Dynamic, 1> U;
+// Matrix<float, Dynamic, 1> 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<int> dims {0, 1, 2};
+ QList<Matrix<float, Dynamic, 1>> results = QtConcurrent::blockingMapped(
+ dims,
+ [&](int j)
+ {
+ Matrix<float, Dynamic, 1> 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++;
}
}
generated by cgit v1.2.3-70-g09d2 (git 2.46.0) at 2025-10-06 04:46:54 +0000