Added final pdfs

BP_Arnaud_Fauconnet.pdf

@@ -0,0 +1 @@
+bachelorproject.pdf

bachelorproject.tex

@@ -5,25 +5,45 @@
 \usepackage[]{subcaption}
 \usepackage[]{float}
 \usepackage[]{multicol}
+\usepackage{mathtools}
+\usepackage[]{listings}
 
 \usepackage{tikz}
 \usepackage{tkz-euclide}
-\usetikzlibrary{external,shapes,through}
+\usetikzlibrary{external,shapes,through,arrows}
 
 \tikzexternalize[prefix=figures/]
 \tikzstyle{none}=[]
 \input{./tikzs/styles.tikzstyles}
 
+\tikzset{>=stealth}
+
 \pgfdeclarelayer{nodelayer}
 \pgfdeclarelayer{edgelayer}
 \pgfsetlayers{edgelayer,nodelayer,main}
 
+
+\lstset{language=C++,
+	basicstyle=\ttfamily,
+	keywordstyle=\color{blue}\ttfamily,
+	stringstyle=\color{red}\ttfamily,
+	commentstyle=\color{green}\ttfamily,
+	morecomment=[l][\color{magenta}]{\#},
+	emph={std, vector, vec2d, polygon, collision},
+	emphstyle=\color{orange}\ttfamily,
+	numbers=left,
+	morekeywords={uint},
+	numberstyle=\footnotesize,
+	numbersep=5pt,
+	frame=lines,
+	% breaklines=true,
+	% breaklines=true,
+	% postbreak=\raisebox{0ex}[0ex][0ex]{\space\ensuremath{\hookrightarrow}},
+}
+
 \graphicspath{{../figures/}{./figures/}}
 
-
 \newcommand*{\vv}[1]{\overrightarrow{#1}}
-
-
 \newcommand*{\figref}[1]{\figurename~\ref{#1}}
 
 \captionsetup{labelfont={bf}}
@@ -47,25 +67,25 @@
 	}
 }
 
-\abstract {
-	Physics engines are a fun and interesting way to learn about a lot of
-	different subjects. First the theoretical concepts, such as the equations
-	that dictate the motion of the objects, together with their components, need
-	to be thoroughly understood. Then there is the necessity of finding a way to
-	represent all of those concepts in a given programming language and to
-	make them as efficient as possible so that the simulation runs fluidly.
-	The task to be completed here was to extend an already existing
-	physics engine that only made circles bounce off each other. The extension
-	was focused on having the ability to generate some arbitrary polygons and
-	make them bounce off each other in a physically accurate way. The main
-	issues that rose up during the development of the extension: determining the
-	inertia of a arbitrary polygon, which is important for realistic
-	impacts; having an accurate collision detection system, which allows the
-	engine to know when to make two polygons bounce off each other. Once those
-	aspects were worked on and polished, the rest of the implementation went
-	smoothly.
-
-}
+\abstract {Physics engines are a fun and interesting way to learn
+  about a the laws of physics, as well as computer science.  They
+  provide a real-time simulation of common physical phenomena, and
+  therefore illustrate theoretical concepts such as the equations that
+  dictate the motion of objects, The goal of this project was to
+  extend an existing physics engine built for demonstration purposes.
+  This engine was initially designed and developed to simulate
+  circular objects (``balls'') in 2D. With this project, we intended
+  to extend this engine to also simulate arbitrary polygons, again in
+  a physically accurate way. The main technical challenges of the
+  project is therefore the correct simulation of the dynamics of
+  rigid, polygonal objects.  In particular, we developed a model of
+  polygonal rigid objects; we implemented a simulation of their
+  inertial motion, possibly in the presence of a constant force field
+  such as gravity; we detect collisions between objects; we compute
+  and then simulate the dynamic effects of collisions.  The
+  simulations are animated and displayed in real-time.  It is also
+  therefore crucial that the simulation code be efficient to obtain
+  smooth animations.}
 
 \counterwithin{figure}{section}
 \counterwithin{equation}{section}

sections/appendix.tex

@@ -3,7 +3,7 @@
 \section{Calculations}
 \label{appendix:calculations}
 
-\paragraph{Moment of inertia of rectangle}
+\subsection{Moment of inertia of rectangle}
 \begin{equation}
 	\label{eq:rect_moment_long}
 	\begin{split}
@@ -22,7 +22,7 @@
 \end{equation}
 
 \newpage
-\paragraph{Moment of inertia of sub-triangle of regular polygon}
+\subsection{Moment of inertia of sub-triangle of regular polygon}
 Before starting the calculations, it is to be noted that according to Figure
 \ref{fig:subtriangle}, we have that
 $$ \tan\left(\frac{\theta}{2}\right) = \frac{\frac{l}{2}}{h} = \frac{l}{2h} $$
@@ -46,7 +46,7 @@ it will be useful to simplify the result of the integral.
 \end{equation}
 
 \newpage
-\paragraph{Moment of inertia of sub-triangle of arbitrary polygon} Recall
+\subsection{Moment of inertia of sub-triangle of arbitrary polygon} Recall
 equation \ref{eq:r} defines
 $$ \vec r = \alpha \vv{CA} + \beta \alpha \vv{AB} $$
 \begin{equation}
@@ -62,3 +62,37 @@ $$ \vec r = \alpha \vv{CA} + \beta \alpha \vv{AB} $$
 		&= \frac{\rho hb}{4} \left(\frac{1}{3}\vv{AB}^2 + \vv{AB} \cdot \vv{CA} + \vv{CA}^2\right) \\
 	\end{split}
 \end{equation}
+
+\newpage
+\subsection{Solving for impulse parameter}
+\label{app:impulse_long}
+We start with equation \ref{eq:vp2n}:
+\begin{equation*}
+	\begin{split}
+		\vec v_{p2} \cdot \vec n &= - e \vec v_{p1} \cdot \vec n\\
+		\left(\vec v_{ap2} - \vec v_{bp2}\right)\cdot \vec n &= - e \vec v_{p1} \cdot \vec n\\
+		\left( \vec v_{a2} + \omega_{a2} \times \vec r_{ap} - \vec v_{b2} - \omega_{b2} \times \vec r_{bp}
+		\right)\cdot \vec n &= - e \vec v_{p1} \cdot \vec n\\
+	\end{split}
+\end{equation*}
+
+We now expand the bracket on the left-hand side using the equations \ref{eq:va2}
+- \ref{eq:omega_b2} and then simplify with equation \ref{eq:vp1}.
+
+\[
+	\begin{split}
+		\left( \vec v_{a1} +  \frac{j\vec n}{m_a} + \omega_{a1} + \frac{\vec r_{ap} \times j\vec n}{I_a} \times \vec r_{ap} - \vec v_{b1} +  \frac{j\vec n}{m_b} - \omega_{b1} + \frac{\vec r_{bp} \times j\vec n}{I_b} \times \vec r_{bp} \right)\cdot \vec n &= - e \vec v_{p1} \cdot \vec n\\
+		\left( \frac{j\vec n}{m_a} + \frac{\vec r_{ap} \times j\vec n}{I_a} \times \vec r_{ap} +  \frac{j\vec n}{m_b} + \frac{\vec r_{bp} \times j\vec n}{I_b} \times \vec r_{bp} \right)\cdot \vec n &= - (1+e) \vec v_{p1} \cdot \vec n\\
+	\end{split}
+\]
+Using the triple scalar product rule, we can derive that
+\[
+	j \left( \frac{1}{m_a} + \frac{1}{m_b} + \frac{\left(\vec r_{ap} \times \vec n\right)^2}{I_a} + \frac{\left( \vec r_{bp} \times \vec n \right)^2}{I_b}
+	\right) = -(1+e) \vec v_{p1} \cdot \vec n
+\]
+and therefore
+\begin{equation}
+	j = \frac{ - (1+e) \cdot \vec v_{ap1} \cdot \vec n }{\frac{1}{m_a} + \frac{1}{m_b} +
+		\frac{\left( \vec r_{ap} \times \vec n \right)^2}{I_a} + \frac{\left( \vec
+			r_{bp} \times \vec n \right)^2}{I_b}}
+\end{equation}

sections/conclusion.tex

@@ -1 +1,26 @@
 \section{Conclusion}
+
+The overall development of this extension was very instructive. Getting to
+discover some theoretical concepts, having to understand them and then being
+able to apply them in the code is a very rewarding experience. The objective of
+project is considered accomplished. In figure \ref{fig:before_after}, we can see
+a comparison of the state of the application before and after the extension.
+
+\begin{figure}[H]
+	\centering
+	\hfill
+	\begin{subfigure}[]{.4\textwidth}
+		\centering
+		\includegraphics[width=.8\textwidth]{before}
+		\caption{Original state of the simulation}
+	\end{subfigure}
+	\hfill
+	\begin{subfigure}[]{.4\textwidth}
+		\centering
+		\includegraphics[width=.8\textwidth]{after}
+		\caption{State after the extension}
+	\end{subfigure}
+	\hfill\null
+	\caption{Before v. After}
+	\label{fig:before_after}
+\end{figure}

sections/implementation.tex

@@ -1 +1,178 @@
 \section{Implementation}
+
+This section will dive into the actual implementing of the project, trying to
+describe how the theoretical concepts in section \ref{sec:theory}.
+
+\subsection{Structure}
+
+The added code to the project can be divided into three parts
+\begin{itemize}
+	\item polygon generator;
+	\item collision detection;
+	\item collision resolution.
+\end{itemize}
+
+Each part will be explained in the following sub-sections.
+
+\subsubsection{Polygon generator}
+Before talking about how we generate polygons, let's first talk about how the
+polygons are defined in this project.
+
+In the file \texttt{polygons.h} we define a polygon as
+
+\begin{lstlisting}[caption={Polygon class (simplified)},label={lst:polygon}]
+class polygon {
+    std::vector<vec2d> points;
+
+    vec2d center;
+    double angle;
+
+    double inertia;
+    double mass;
+
+    std::vector<vec2d> global_points();
+
+    vec2d speed;
+    double angular_speed;
+}
+\end{lstlisting}
+
+Polygons possesses multiple fields. Firstly we have
+\lstinline{std::vector<vec2d> points}, a collection of
+\lstinline{vec2d} objects, which represent the set of ordered points that
+compose the polygon. Those points are expressed in local coordinates, and the
+center of mass of the polygon is placed the origin.
+
+Now that we know how the polygon is composed, we want to move it around the
+space it lives in, that is the purpose of the \lstinline{vec2d center} field.
+It represents where the center of mass of the polygon is located in simulation.
+Since the shapes also rotate, we need to keep track of the rotation of the
+polygon, hence the use of \lstinline{double angle}. All of these fields are used
+when computing the result of the method \lstinline{std::vector<vec2d> global_points()}.
+
+Finally, we have \lstinline{double inertia} and \lstinline{double mass} which,
+as their name suggest, store the value for the polygon's inertia and mass, that
+are used to calculate the polygons final speed and angular speed, who are
+represented by \lstinline{vec2d speed} and \lstinline{double angular_speed}.
+
+\paragraph{Polygon generator} Now that we know how polygons are represented in
+our simulation, we can generate them. In \texttt{polygon\_generator.h} (and its
+implementation file \texttt{polygon\_generator.cc}), there are some functions
+for that. In Listing \ref{lst:polygon_gen} we can see the signature of the
+functions that generate
+\begin{itemize}
+	\item a rectangle of a certain width and height;
+	\item a square of a certain side length;
+	\item a triangle given two side lengths and an angle between them;
+	\item a regular polygon;
+	\item an arbitrary polygon.
+\end{itemize}
+
+\begin{lstlisting}[caption={Polygon Generator header
+file},label={lst:polygon_gen}]
+namespace poly_generate {
+    polygon rectangle(double width, double height);
+
+    inline polygon square(double width) {
+        assert(width > 0);
+        return rectangle(width, width, label);
+    };
+
+    polygon triangle(double side1, double side2, double angle);
+
+    polygon regular(double radius, uint n_sides);
+
+    polygon general(std::vector<vec2d> points);
+};
+\end{lstlisting}
+
+The implementation of those functions are fairly straight forward, they generate
+the appropriate number of points in the correct place. After what they calculate
+the mass and subsequent inertia of the polygon as shown in section
+\ref{sub:moment}.
+
+\subsubsection{Collision}
+\label{sub:implementation-collision-detection}
+
+The algorithm described in section~\ref{sub:vertex-collision} is implemented in
+\texttt{collision.cc} and exposed to other modules with \texttt{collision.h}.
+
+The module exposes a collision structure and a collides function.
+
+\begin{lstlisting}[caption={Collision header file},label={lst:collision}]
+struct collision {
+    bool collides = false;
+    vec2d impact_point;
+    vec2d n;
+    vec2d overlap;
+};
+
+extern collision collides(polygon& p, polygon& q);
+\end{lstlisting}
+
+The \lstinline{collides} function takes in two polygons and checks with
+vertex-collision algorithm whether they collide or not. The result is return
+through an instance of \lstinline{struct collision}. The \lstinline{bool collides}
+and \lstinline{vec2d impact_point} fields is self-explanatory. The
+\lstinline{vec2d n} is the normal vector that pushes \lstinline{polygon& p} away
+from \lstinline{polygon& q}. Finally, \lstinline{vec2d overlap} is a scalar
+multiplication of \lstinline{vec2d n}. Where the latter is a normalized vector,
+the former represents how deep \lstinline{p} is in \lstinline{q}. If we push
+\lstinline{p} the exact amount of \lstinline{overlap}, then \lstinline{p} would
+just be touching \lstinline{q} with point \lstinline{impact_point}, and no
+further overlap would occur.
+
+In reality, we do not simply push \lstinline{p} by \lstinline{overlap}, but we
+push \lstinline{p} and \lstinline{q} away from each other proportionally to
+their mass.
+
+\paragraph{Collision resolution} In \texttt{polygons.cc}, among other functions,
+we apply the collision resolution, which simply uses the physics results found
+in \ref{sub:resolution}.
+
+\subsection{Optimization}
+
+In order to optimise the collision detection and to avoid having to do a check
+whether each vertex of every polygon was colliding with each edge of every other
+polygon we decided to apply the bounding box acceleration structure.
+
+The bounding box of an object is, in two dimensions, the rectangle that contains
+the object and which the sides are aligned with axis of the coordinate system.
+
+\begin{figure}[H]
+	\centering
+	\inputtikz[.7]{bounding_box}
+	\caption{Bounding box example}
+	\label{fig:bounding-box}
+\end{figure}
+
+To get the four coordinates that compose, we just perform a linear scan through
+all the points that compose a polygon, record the minimum and maximum of both
+they $x$ and $y$ coordinate.
+
+Once those points have been found, we just perform some basic axis-aligned
+rectangle collision detection, which is much less computationally expensive than
+checking each vertex-edge pair.
+
+The complexity of the collision detection algorithm went from $\mathcal{O}(n^2)$
+to $\mathcal{O}(n)$, where $n$ is the total number of vertices across all
+polygons in the simulation.
+
+
+\subsection{Known issues}
+
+The simulation has one major flaw, and it is inherent to the last part of
+section \ref{sub:implementation-collision-detection} where we talk about the
+\lstinline{overlap} vector: we said that this vector allows to displace both
+polygons so that and the end of the collision resolution, they are not
+overlapping anymore. This was done to avoid issues where the collision
+resolution had taken place in one frame, but the polygons were still overlapping
+in the following frame, which lead to polygons getting stuck together.
+
+The unfortunate consequence of such a decision is that, when we start playing
+around with the restitution coefficient and all the shapes start falling to the
+ground, the shapes are not able to rest still. Since there is gravity, at each
+frame their speed get updated to simulate its effect, but since they are lying
+on the floor, the collision detection algorithm kicks in right away, resulting
+in the polygon getting pushed completely up, since the ground has infinite mass.
+This results in the polygon bouncing on the floor instead of lying down.

sections/theoretical_background.tex

@@ -1,44 +1,61 @@
 \section{Theoretical Background}
+\label{sec:theory}
 The theoretical background is everything related to the physics part of the
-project. It covers the calculating the inertia of different types of polygons;
-different algorithms to detect whether there is a collision between two
-polygons; the resolution of the collision, i.e. finding the final
+project. It covers the calculation of the inertia of different types of polygons;
+different algorithms to detect whether there is a collision between
+any two polygons; and the resolution of the collision, that is, finding the final
 velocity vectors and angular speed of those polygons.
 
-\subsection{Moment of inertia}
+\subsection{Moment of Inertia}
+\label{sub:moment}
 
-The inertia of an object refers to the tendency of an object to resist a change
-of its state of motion or rest, it describes how the object behaves when forces
-are applied to it. An object with a lot of inertia requires more force to change
-its motion, either to make it move if it's at rest or to stop it if it's already
-moving. On the other hand, an object with less inertia is easier to set in
-motion or bring to a halt.
+The inertia of an object refers to the tendency of an object to resist
+a change of its state of motion or rest, it describes how the object
+behaves when forces are applied to it. In particular, the mass of a
+point object (inertial mass) determines the inertia of that objects.
+Thus an object with a large mass requires more force to change its
+motion, either to make it move from an initial state at rest or to
+stop it if it is already moving. On the other hand, an object with a
+smaller mass is easier to set in motion or to bring to a halt.
 
-The moment of inertia is similar but is used in a slightly different context, it
-specifically refers to the rotational inertia of an object. It measures an
-object's resistance to changes in its rotational motion and how its mass is
-distributed with respect to is axis of rotation.
+For objects that are not point-like, the motion includes a rotational
+component.
 
-In the case of this project the axis of rotation is the one along the $z$-axis
-(perpendicular to the plane of the simulation) and placed at the barycenter of
-the polygon.
+The moment of inertia is analogous to the mass in determining the
+inertia of an object with respect to its rotational motion. In other
+words, the moment of inertia determines the rotational inertia of an
+object, that is, the object's resistance to changes in its rotational
+motion.  The moment of inertia of an object is therefore computed with
+respect to a rotation axis, and it is determined by the way the mass
+of the object is distributed with respect to the chosen axis of
+rotation.
 
-The general formula for the moment of inertia is
+
+For the purpose of this project, we are interested in a 2D simulation.
+This means that the axis of rotation is always parallel to the
+$z$-axis, perpendicular to the plane of the simulation.  Furthermore,
+we model the motion of an object using the barycenter of an object as
+its reference point.  We therefore compute the moment of inertia of a
+polygon with respect to a rotation axis that is perpendicular to the
+2D plane and that goes through the center of mass of the polygon.
+
+In general, the moment of inertia about a rotation axis $Z$ of a point
+mass $m$ at distance $r$ from $Z$ is $I_m=mr^2$.  Therefore, for a
+generic 2D object $Q$, we integrate over the whole 2D object:
 \begin{equation}
 	\label{eq:moment_general}
-	I_Q = \int \vec r^2 \rho(\vec r) \diff \mathcal{A}
+	I_Q = \int_{\mathcal{A}}\|\vec r\|^2 \rho(\vec r) \diff \mathcal{A}
 \end{equation}
-where $\rho$ is the density of object $Q$ in the point $\vec r$ across the
-small pieces of area $\mathcal A$ of the object.
+where $\rho$ is the density of object $Q$ (mass per unit area) in the
+point $\vec r$ from the across the small pieces of area $\mathcal A$
+of the object.
 
-In our case, since we are implementing a 2D engine we can use the $\mathbb{R}^2$
-coordinate systems, thus the formula becomes
-$$ I_Q = \iint \rho(x, y) \vec r^2 \diff x\diff y$$
-and since the requirements express that the mass of the polygons is spread
-uniformly across its surface, the formula finally becomes
+In our case, we use 2D Cartesian coordinates centered in the point
+(axis) about which we compute he moment of inertia.  Furthermore, we
+consider polygons of uniform density.  We can therefore compute:
 \begin{equation}
 	\label{eq:moment}
-	I_Q = \rho \iint x^2 + y^2 \diff x\diff y
+	I_Q = \rho \iint(x^2 + y^2)\diff x\diff y
 \end{equation}
 
 The bounds of the integral depend on the shape of the polygon. In the following
@@ -264,6 +281,7 @@ to the overall polygon.
 where, $P_{n+1} = P_1$ in the case of $i = n$.
 
 \subsection{Collision detection}
+\label{sub:collision-detection}
 
 Collision detection, as the name suggests, are the algorithms used to detect
 whether two polygons are colliding. The result of this procedure must be an
@@ -324,6 +342,7 @@ algorithm of our own. Moreover, SAT only supports convex polygons, which limits
 the original objective of the project, which was to have any arbitrary polygon.
 
 \subsubsection{Vertex collisions}
+\label{sub:vertex-collision}
 
 The solution that was adopted for the project, after trying SAT, was a more
 intuitive one, developed by Prof. Carzaniga. The idea is simple: check if a
@@ -462,3 +481,164 @@ find the normal. The results look realistic enough to be accepted.
 
 \subsection{Collision resolution}
 \label{sub:resolution}
+
+The collision resolution is the last step in the processing of the collision.
+Algorithmically, it is much less heavy than collision detection, since, once the
+simulation has two colliding polygons, a point of impact and a normal vector,
+it's just a case of applying the rigid body physics formulas to the polygons
+that are colliding. This part has been helped a lot by the works of Erik Neumann
+\cite{collision:resolution-site} and Chris Hecker
+\cite{collision:resolution-paper}.
+
+\begin{figure}[H]
+	\centering
+	\inputtikz[.5]{collision_resolution}
+	\caption{Collision resolution between polygons $A$ and $B$}
+	\label{fig:collision_resolution}
+\end{figure}
+
+\paragraph{Variable definition} Before getting into any maths, let's define some
+variables that we are going to use
+
+\begin{multicols}{2}
+	\begin{itemize}
+		\item $m_a, m_b = $ mass of the bodies $A$ and $B$
+		\item $\vec r_{ap} = $ distance vector from center of mass of body $A$ to
+		      point $P$
+		\item $\vec r_{bp} = $ distance vector from center of mass of body $A$ to
+		      point $P$
+		\item $\omega_{a1}, \omega_{b1} = $ initial angular velocity of bodies
+		      $A, B$
+		\item $\omega_{a2}, \omega_{b2} = $ final angular velocity of bodies
+		      $A, B$
+		\item $\vec v_{a1}, \vec v_{b1} =$ initial velocities of center of mass
+		      bodies $A, B$
+		\item $\vec v_{a2}, \vec v_{b2} =$ final velocities of center of mass
+		      bodies $A, B$
+		\item $\vec v_{ap1}=$ initial velocity of impact point $P$ on body $A$
+		\item $\vec v_{bp1}=$ initial velocity of impact point $P$ on body $B$
+		\item $\vec v_{p1}=$ initial relative velocity of impact points on body $A, B$
+		\item $\vec v_{p2}=$ final relative velocity of impact points on body $A, B$
+		\item $\vec n=$ normal vector
+		\item $ e=$ elastic coefficient (0 = inelastic, 1 = perfectly elastic)
+	\end{itemize}
+\end{multicols}
+
+Some of those variables, like the ones in the left column, are already given by
+the simulation, some of them have been computed thanks to the algorithms in
+section \ref{sub:collision-detection} (normal vector and position of point $P$),
+and some have still to be defined mathematically, such as $\vec v_{ap1},\vec
+	v_{bp1},\vec v_{p1}$ and $\vec v_{p2}$. So the velocity vectors of impact
+point $P$ on both bodies, before the collision, are
+\begin{equation} \label{eq:vabp1}
+	\begin{split}
+		\vec v_{ap1} = \vec v_{a1} + \omega_{a1} \times \vec r_{ap} \\
+		\vec v_{bp1} = \vec v_{b1} + \omega_{b1} \times \vec r_{bp}
+	\end{split}
+\end{equation}
+
+Similarly, the final velocities are
+\begin{equation} \label{eq:vabp2}
+	\begin{split}
+		\vec v_{ap2} = \vec v_{a2} + \omega_{a2} \times \vec r_{ap}\\
+		\vec v_{bp2} = \vec v_{b2} + \omega_{b2} \times \vec r_{bp}
+	\end{split}
+\end{equation}
+Here we are regarding the angular velocity as a 3-dimensional vector
+perpendicular to the plane, so that the cross product is calculated as
+
+$$ \omega \times \vec r = \begin{pmatrix} 0\\0\\\omega \end{pmatrix} \times
+	\begin{pmatrix} r_x\\r_y\\0 \end{pmatrix} = \begin{pmatrix} -\omega r_y \\
+		\omega r_x  \\0\end{pmatrix}  $$
+
+We these variables, we can finally define the relative velocities $\vec
+	v_{p1}$ and $\vec v_{p2}$
+
+\begin{equation}
+	\label{eq:vp1}
+	\begin{split}
+		\vec v_{p1} = \vec v_{ap1} - \vec v_{bp2}\\
+		\vec v_{p2} = \vec v_{ap2} - \vec v_{bp2}
+	\end{split}
+\end{equation}
+If we expand by using \ref{eq:vabp1} and \ref{eq:vabp2}, we get
+\begin{equation}
+	\begin{split}
+		\vec v_{p1} = \vec v_{a1} + \omega_{a1} \times \vec r_{ap} - \vec v_{b1}
+		+ \omega_{b1} \times \vec r_{bp}\\
+		\vec v_{p2} = \vec v_{a2} + \omega_{a2} \times \vec r_{ap} - \vec v_{b2} + \omega_{b2} \times \vec r_{bp}
+	\end{split}
+\end{equation}
+
+The relative velocity of point $P$ along the normal vector $\vec n$ is
+$$ \vec v_{p1} \cdot \vec n $$
+Note that for a collision to occur this relative normal velocity must be
+negative (that is, the objects must be approaching each other). Let e be the
+elasticity of the collision, having a value between 0 (inelastic) and 1
+(perfectly elastic). We now make an important assumption in the form of the
+following relation
+\begin{equation}
+	\label{eq:vp2n}
+	\vec v_{p2} \cdot \vec n = - e \vec v_{p1} \cdot \vec n
+\end{equation}
+This says that the velocity at which the objects fly apart is proportional to
+the velocity with which they were coming together. The proportionality factor is
+the elasticity $e$.
+
+\paragraph{Collision Impulse} In simple terms, collision impulse refers to the
+change in momentum experienced by an object during a collision. It is a measure
+of the force applied to an object over a short period of time. We imagine that
+during the collision there is a very large force acting for a very brief period
+of time. If you integrate (sum) that force over that brief time, you get the
+impulse.
+
+In our simulation, we are assuming that no friction is happening during the
+collision, so that the only force we have to consider is the one of along the
+normal vector $\vec n$. The friction would create a force perpendicular to the
+normal, and it would make things a little too complicated for the scope of this
+project.
+
+Since the only force we consider is the normal one, we can consider the net
+impulse to be $j \vec n$, where $j$ is the impulse parameter. The body $A$ will
+experience the net impulse $j \vec n$ and body $B$ will experience it's negation
+$- j \vec n$ since the force that $B$ experience is equal an opposite to the one
+experienced by $A$. The impulse is a change in momentum. Momentum has units of
+velocity times mass, so if we divide the impulse by the mass we get the change
+in velocity. We can relate initial and final velocities as
+
+\begin{equation}
+	\label{eq:va2}
+	\vec v_{a2} = \vec v_{a1} +  \frac{j\vec n}{m_a}
+\end{equation}
+\begin{equation}
+	\label{eq:vb2}
+	\vec v_{b2} = \vec v_{b1} -  \frac{j\vec n}{m_b}
+\end{equation}
+
+
+For the final angular speed, it's change from the impulse $j\vec n$ is given by
+$\vec r_{ap} \times j \vec n$. We can divide the result by the moment of
+inertia, which was calculated in section \ref{sub:moment}, to get convert the
+change in angular momentum into change in angular speed. Similarly as above, we
+can relate the initial and final angular velocity as
+
+\begin{equation}
+	\label{eq:omega_a2}
+	\omega_{a2} = \omega_{a1} + \frac{\vec r_{ap} \times j\vec n}{I_a}
+\end{equation}
+\begin{equation}
+	\label{eq:omega_b2}
+	\omega_{b2} = \omega_{b1} - \frac{\vec r_{bp} \times j\vec n}{I_b}
+\end{equation}
+
+\paragraph{Solving for the impulse parameter}
+Now we have everything we need to solve for the impulse parameter $j$. The bulk
+of the calculations can be found in section \ref{app:impulse_long} of the
+appendix \ref{appendix:calculations}. But basically we start with equation
+\ref{eq:vp2n} and we end up with
+\begin{equation}
+	\label{eq:j}
+	j = \frac{ - (1+e) \cdot \vec v_{ap1} \cdot \vec n }{\frac{1}{m_a} + \frac{1}{m_b} +
+		\frac{\left( \vec r_{ap} \times \vec n \right)^2}{I_a} + \frac{\left( \vec
+			r_{bp} \times \vec n \right)^2}{I_b}}
+\end{equation}

tikzs/arbitrary_subtriangle.tikz

@@ -37,7 +37,7 @@
 		\draw [style=none] (15.center) to (14.center);
 		\draw [style=none] (14.center) to (13.center);
 		\draw [style=none] (13.center) to (12.center);
-		\draw [style=Axis] (0.center) to (16);
+		\draw [style=Vector] (0.center) to (16);
 		\draw [style=Dotted] (0.center) to (21);
 		\draw (23.center) to (26.center);
 		\draw (26.center) to (25.center);

tikzs/bounding_box.tikz

@@ -0,0 +1,71 @@
+\begin{tikzpicture}
+	\begin{pgfonlayer}{nodelayer}
+		\node [style=none] (0) at (0, -1) {};
+		\node [style=none] (1) at (0, 8) {};
+		\node [style=none] (2) at (-1, 0) {};
+		\node [style=none] (3) at (15, 0) {};
+		\node [style=none] (4) at (4, 4) {};
+		\node [style=none] (5) at (2, 1) {};
+		\node [style=none] (6) at (6, 3) {};
+		\node [style=none] (7) at (9, 3) {};
+		\node [style=none] (8) at (11, 2) {};
+		\node [style=none] (9) at (10.5, 6.5) {};
+		\node [style=none] (10) at (12.75, 5.75) {};
+		\node [style=none] (11) at (13.5, 3.5) {};
+		\node [style=none] (12) at (8, 5) {};
+		\node [style=none] (13) at (2, 4) {};
+		\node [style=none] (14) at (6, 1) {};
+		\node [style=none] (15) at (6, 4) {};
+		\node [style=none] (16) at (8, 2) {};
+		\node [style=none] (17) at (8, 6.5) {};
+		\node [style=none] (18) at (13.5, 6.5) {};
+		\node [style=none] (19) at (13.5, 2) {};
+		\node [style=none] (20) at (2, 0) {};
+		\node [style=none] (21) at (6, 0) {};
+		\node [style=none] (22) at (0, 4) {};
+		\node [style=none] (23) at (0, 1) {};
+		\node [style=none] (24) at (-0.25, 1) {$A_{y_{\min}}$};
+		\node [style=none] (25) at (-0.25, 4) {$A_{y_{\max}}$};
+		\node [style=none] (26) at (2, -0.25) {$A_{x_{\min}}$};
+		\node [style=none] (27) at (6, -0.25) {$A_{x_{\max}}$};
+		\node [style=none] (28) at (8, 0) {};
+		\node [style=none] (29) at (8, -0.25) {$B_{x_{\min}}$};
+		\node [style=none] (30) at (13.5, 0) {};
+		\node [style=none] (31) at (13.5, -0.25) {$B_{x_{\max}}$};
+		\node [style=none] (32) at (0, 2) {};
+		\node [style=none] (33) at (-0.25, 2) {$B_{y_{\min}}$};
+		\node [style=none] (34) at (0, 6.5) {};
+		\node [style=none] (35) at (-0.25, 6.5) {$B_{y_{\max}}$};
+		\node [style=none] (36) at (4, 2.75) {$A$};
+		\node [style=none] (37) at (10.5, 4.5) {$B$};
+	\end{pgfonlayer}
+	\begin{pgfonlayer}{edgelayer}
+		\draw [style=Vector] (0.center) to (1.center);
+		\draw [style=Vector] (2.center) to (3.center);
+		\draw (5.center) to (4.center);
+		\draw (4.center) to (6.center);
+		\draw (6.center) to (5.center);
+		\draw (12.center) to (9.center);
+		\draw (9.center) to (10.center);
+		\draw (10.center) to (11.center);
+		\draw (11.center) to (8.center);
+		\draw (8.center) to (7.center);
+		\draw (7.center) to (12.center);
+		\draw [style=Dashed] (13.center) to (5.center);
+		\draw [style=Dashed] (5.center) to (14.center);
+		\draw [style=Dashed] (14.center) to (15.center);
+		\draw [style=Dashed] (15.center) to (13.center);
+		\draw [style=Dashed] (16.center) to (17.center);
+		\draw [style=Dashed] (17.center) to (18.center);
+		\draw [style=Dashed] (18.center) to (19.center);
+		\draw [style=Dashed] (19.center) to (16.center);
+		\draw [style=Dotted] (19.center) to (30.center);
+		\draw [style=Dotted] (16.center) to (28.center);
+		\draw [style=Dotted] (16.center) to (32.center);
+		\draw [style=Dotted] (13.center) to (22.center);
+		\draw [style=Dotted] (17.center) to (34.center);
+		\draw [style=Dotted] (5.center) to (20.center);
+		\draw [style=Dotted] (5.center) to (23.center);
+		\draw [style=Dotted] (14.center) to (21.center);
+	\end{pgfonlayer}
+\end{tikzpicture}

tikzs/isosceles.tikz

@@ -17,8 +17,8 @@
 		\node [style=none] (14) at (2, .5) {\Large$\frac{\theta}{2}$};
 	\end{pgfonlayer}
 	\begin{pgfonlayer}{edgelayer}
-		\draw [style=Axis] (0.center) to (1.center);
-		\draw [style=Axis] (2.center) to (3.center);
+		\draw [style=Vector] (0.center) to (1.center);
+		\draw [style=Vector] (2.center) to (3.center);
 		\draw (4.center) to (5.center);
 		\draw (5.center) to (6.center);
 		\draw (6.center) to (4.center);

tikzs/rectangle_inertia2d.tikz

@@ -25,7 +25,7 @@
 		\draw (1.center) to (2.center);
 		\draw (2.center) to (3.center);
 		\draw (3.center) to (0.center);
-		\draw [style=Axis] (16.center) to (11.center);
-		\draw [style=Axis] (15.center) to (12.center);
+		\draw [style=Vector] (16.center) to (11.center);
+		\draw [style=Vector] (15.center) to (12.center);
 	\end{pgfonlayer}
 \end{tikzpicture}

tikzs/rectangle_inertia3d.tikz

@@ -25,8 +25,8 @@
 		\draw (1.center) to (2.center);
 		\draw (2.center) to (3.center);
 		\draw (3.center) to (0.center);
-		\draw [style=Axis] (16.center) to (11.center);
-		\draw [style=Axis] (15.center) to (12.center);
-		\draw [style=Axis] (21.center) to (22.center);
+		\draw [style=Vector] (16.center) to (11.center);
+		\draw [style=Vector] (15.center) to (12.center);
+		\draw [style=Vector] (21.center) to (22.center);
 	\end{pgfonlayer}
 \end{tikzpicture}