add videos, concurrency & add dirty init

11 files changed, 244 insertions, 154 deletions

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
--- a/readme-videos/armadillo.gif
+++ b/readme-videos/armadillo.gif
diff --git a/readme-videos/bean.gif b/readme-videos/bean.gif
index 5a80c0b..be92079 100644
--- a/readme-videos/bean.gif
+++ b/readme-videos/bean.gif
diff --git a/readme-videos/bunny-conc.gif b/readme-videos/bunny-conc.gif
new file mode 100644
index 0000000..459c927
--- /dev/null
+++ b/readme-videos/bunny-conc.gif
diff --git a/readme-videos/bunny-norm.gif b/readme-videos/bunny-norm.gif
new file mode 100644
index 0000000..4774bbb
--- /dev/null
+++ b/readme-videos/bunny-norm.gif
diff --git a/readme-videos/peter.gif b/readme-videos/peter.gif
index bb00f4b..daa6ce2 100644
--- a/readme-videos/peter.gif
+++ b/readme-videos/peter.gif
diff --git a/readme-videos/sphere.gif b/readme-videos/sphere.gif
index 47a131a..e0a8ae1 100644
--- a/readme-videos/sphere.gif
+++ b/readme-videos/sphere.gif
diff --git a/readme-videos/teapot.gif b/readme-videos/teapot.gif
index 7c202bb..e89fcba 100644
--- a/readme-videos/teapot.gif
+++ b/readme-videos/teapot.gif
diff --git a/readme-videos/tetrahedron.gif b/readme-videos/tetrahedron.gif
index 866f1e2..0a8f139 100644
--- a/readme-videos/tetrahedron.gif
+++ b/readme-videos/tetrahedron.gif
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++;
 				}
 			}
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 <Eigen/Core>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
+#include <QList>
+#include <QtConcurrent>
 
 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<float, 3, Eigen::Dynamic> PM; // position matrix
+	typedef Eigen::Matrix<float, 3, 3> RM; // rotation matrix
+	typedef Eigen::SparseMatrix<float, Eigen::RowMajor> Sparse;
+
+	std::vector<RM> m_rotations;
+	Sparse m_edge_weights;
+	PM m_b, m_b_fixed;
+	std::vector<Eigen::Index> m_vtx_to_free_vtx;
+	Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>> m_solver;
+
+	int m_num_iterations;
+	const char * m_mesh_path;
 };