improve performance and udpate some videos

5 files changed, 184 insertions, 194 deletions

diff --git a/README.md b/README.md
index c637b9a..8642ed7 100644
--- a/README.md
+++ b/README.md
@@ -90,7 +90,7 @@ For this project, we ask that you demonstrate to us that your program achieves t
 ### 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.
+See line 286 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.
 
diff --git a/readme-videos/bean.gif b/readme-videos/bean.gif
index be92079..968d353 100644
--- a/readme-videos/bean.gif
+++ b/readme-videos/bean.gif
diff --git a/readme-videos/peter.gif b/readme-videos/peter.gif
index daa6ce2..161f78a 100644
--- a/readme-videos/peter.gif
+++ b/readme-videos/peter.gif
diff --git a/src/arap.cpp b/src/arap.cpp
index f4ac929..0dcaefd 100644
--- a/src/arap.cpp
+++ b/src/arap.cpp
@@ -6,8 +6,9 @@
 #include <map>
 #include <vector>
 
-#include <QList>
-#include <QtConcurrent>
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
 
 using namespace std;
 using namespace Eigen;
@@ -16,49 +17,43 @@ 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(m_mesh_path, vertices, triangles)) {
-        m_shape.init(vertices, triangles);
-    }
-
-    // Students, please don't touch this code: get min and max for viewport stuff
-    MatrixX3f all_vertices = MatrixX3f(vertices.size(), 3);
-    int i = 0;
-    for (unsigned long i = 0; i < vertices.size(); ++i) {
-        all_vertices.row(i) = vertices[i];
-    }
-    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());
+	m_mesh_path = "meshes/peter.obj";
+
+	vector<Vector3f> vertices;
+	vector<Vector3i> triangles;
+
+	// If this doesn't work for you, remember to change your working directory
+	if (MeshLoader::loadTriMesh(m_mesh_path, vertices, triangles)) {
+		m_shape.init(vertices, triangles);
+	}
+
+	// Students, please don't touch this code: get min and max for viewport stuff
+	MatrixX3f all_vertices = MatrixX3f(vertices.size(), 3);
+	int i = 0;
+	for (unsigned long i = 0; i < vertices.size(); ++i) {
+		all_vertices.row(i) = vertices[i];
+	}
+	coeffMin = all_vertices.colwise().minCoeff();
+	coeffMax = all_vertices.colwise().maxCoeff();
 }
 
+
+
 // Move an anchored vertex, defined by its index, to targetPosition
 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::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;
 
-	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)
 	{
@@ -71,158 +66,153 @@ void ARAP::move(int vertex, Vector3f targetPosition)
 		p_prime.col(i) = new_vertices[i];
 	}
 
-	if (new_num_anchors != num_anchors)
+	// 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)
 	{
-		// 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 (!anchors.contains(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++;
-				}
-			}
+			// 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 (m_vtx_to_free_vtx[i] == -1)
-				{
-					// if we're an anchor, set this position
-					p_prime.col(i) = new_vertices[i];
-				}
-			}
+	// 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];
+		}
+	}
 
-			// 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();
+	// 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();
 
-			// setup bfixed
-			m_b_fixed = PM::Zero(3, count_free);
-			m_b_fixed.resize(3, count_free);
+	// setup bfixed
+	PM b_fixed = PM::Zero(3, count_free);
+	// b_fixed.resize(3, count_free);
 
-			// 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++)
-			{
-				Index free_index = m_vtx_to_free_vtx[i];
-				if (free_index == -1)
-				{
-					// do not add constraints to the linear system
-					continue;
-				}
+	// 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;
+		}
 
-				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));
-				}
+		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)
+			{
+				// if the neighbor is an anchor, add to b_fixed
+				Vector3f location = new_vertices[j];
+				b_fixed.col(free_index) += wij * location;
 			}
-			// create and solve the L matrix
-			L.setFromTriplets(triplets_L.begin(), triplets_L.end());
-			m_solver.compute(L);
+			else
+			{
+				triplets_L.push_back(Tri(free_index, free_index_j, -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());
+	SimplicialLDLT<SparseMatrix<float>> solver;
+	solver.compute(L);
 
-	// do the iterations
 	for (int iter = 0; iter < m_num_iterations; iter++)
 	{
 		// estimate the rotations
-		m_rotations.clear();
-		m_rotations.resize(m_shape.getVertices().size(), RM::Identity(3, 3));
-
+		rotations.clear();
+		rotations.resize(m_shape.getVertices().size(), R::Identity(3, 3));
 		typedef Matrix<float, 3, 3> S;
-		for (Index i = 0; i < m_edge_weights.outerSize(); ++i)
+		for (Index i = 0; i < 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(m_edge_weights, i); it; ++it)
+			for (SparseMatrix<float, RowMajor>::InnerIterator it(edge_weights, i); it; ++it)
 			{
 				// get the weight
 				float wij = it.value();
@@ -243,35 +233,35 @@ void ARAP::move(int vertex, Vector3f targetPosition)
 
 			S I = S::Identity(3, 3);
 			I(2, 2) = (v * u_trans).determinant();
-			m_rotations[i] = v * I * u_trans;
+			rotations[i] = v * I * u_trans;
 		}
 
 		// estimate the positions
 		// make a copy of b_fixed to modify
-		m_b = m_b_fixed;
-		for (Index i = 0; i < m_edge_weights.outerSize(); ++i)
+		PM b = b_fixed;
+		for (Index i = 0; i < edge_weights.outerSize(); ++i)
 		{
 			// TODO add constraint check here
-			Index free_index = m_vtx_to_free_vtx[i];
+			Index free_index = vtx_to_free_vtx[i];
 			if (free_index == -1)
 			{
 				// do not consider constraints
 				continue;
 			}
 
-			for (SparseMatrix<float, RowMajor>::InnerIterator it(m_edge_weights, i); it; ++it)
+			for (SparseMatrix<float, RowMajor>::InnerIterator it(edge_weights, i); it; ++it)
 			{
 				// get the weight
 				float wij = it.value();
 
 				// get the rotation matrices for the vertices
-				const S &R = m_rotations[i];
-				const S &Rj = m_rotations[it.col()];
+				const S &R = rotations[i];
+				const S &Rj = 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;
-				m_b.col(free_index) += pt;
+				b.col(free_index) += pt;
 			}
 		}
 
@@ -292,14 +282,14 @@ void ARAP::move(int vertex, Vector3f targetPosition)
 //			}
 //		}
 
-//		// EXTRA FEATURE - parallelize the solving of the linear system by dimension
+		// 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());
+				U = solver.solve(b.row(j).transpose());
 				return U;
 			});
 
@@ -309,7 +299,7 @@ void ARAP::move(int vertex, Vector3f targetPosition)
 			Index idx = 0;
 			for (int k = 0; k < m_shape.getVertices().size(); ++k)
 			{
-				if (m_vtx_to_free_vtx[k] != -1)
+				if (vtx_to_free_vtx[k] != -1)
 				{
 					p_prime(j, k) = results[j](idx);
 					idx++;
@@ -325,18 +315,18 @@ void ARAP::move(int vertex, Vector3f targetPosition)
 		new_vertices[i] = p_prime.col(i);
 	}
 
-    // Here are some helpful controls for the application
-    //
-    // - You start in first-person camera mode
-    //   - WASD to move, left-click and drag to rotate
-    //   - R and F to move vertically up and down
-    //
-    // - C to change to orbit camera mode
-    //
-    // - Right-click (and, optionally, drag) to anchor/unanchor points
-    //   - Left-click an anchored point to move it around
-    //
-    // - Minus and equal keys (click repeatedly) to change the size of the vertices
-
-    m_shape.setVertices(new_vertices);
-}
+	// Here are some helpful controls for the application
+	//
+	// - You start in first-person camera mode
+	//   - WASD to move, left-click and drag to rotate
+	//   - R and F to move vertically up and down
+	//
+	// - C to change to orbit camera mode
+	//
+	// - Right-click (and, optionally, drag) to anchor/unanchor points
+	//   - Left-click an anchored point to move it around
+	//
+	// - Minus and equal keys (click repeatedly) to change the size of the vertices
+
+	m_shape.setVertices(new_vertices);
+}
\ No newline at end of file
diff --git a/src/arap.h b/src/arap.h
index 10a2904..82eb10e 100644
--- a/src/arap.h
+++ b/src/arap.h
@@ -65,4 +65,4 @@ public:
 
 	int m_num_iterations;
 	const char * m_mesh_path;
-};
+};
\ No newline at end of file