Category Archives: 2D

Printing Pretty Ascii Trees


I’m terrible at writing recursive algorithms. Absolutely horrible. So for practice I decided to try printing a binary tree where all branches have two children, all in the console. It turns out this is a pretty useful quick-n-dirty debugging tool when creating trees.

I really like formatting of Tree for linux so I decided to steal it.

After 3 straight hours of trying out stupid stuff I finally come to a working solution (my earlier attempts were around 150 lines of code, the solution was like 70):

Which will print cool trees like this one:

I then realized that I could use singly linked lists to represent trees with any number of children, which is totally necessary for the parser I am currently writing. So after 30 minutes I was able to form version 2 for arbitrary trees; this version is surprisingly right around the same length and complexity as the previous version:

And will print trees like this monster:

I also ended up finding a pretty cool implementation by a random stack-overflow user that prints trees in horizontal format (which turns out to be much more limited for the purposes of a console, since it can’t print larger trees, but might look a lot prettier for small trees):

It may be preferable to upgrade to the limited edition Box Drawing characters for best viewing pleasure:



In which case you’ll probably want to use some specific character codes (I used 195, 196 and 192 in the above image).


Interactive Cubic Bezier Splines


Cubic Bezier splines look rather beautiful. Splines in general look great and are super useful for all sorts of cool stuff.

Creating splines usually involves defining some control points and moving those around until the desired effect is achieved. However, clicking on or interacting with the spline itself turns out to be a lot harder than just messing with control points.

In particular I want to talk about cubic Bezier splines. This kind of spline is defined as 3 * n + 1 control points, where n is the number of Bezier curves in the spline. The first four points represent the first curve, then the last point and the next three represent the next curve, and so on.

The resulting spline can look like the curved arrows in this graph:

Actually interacting with these splines would be pretty difficult since cartesian coordinates are not in the same basis as the spline itself. What is needed is the ability to express a relation between these two spaces.

Since I haven’t studied curves in great detail there’s much math that goes beyond me (for the time being) in solving this problem. However in Graphics Gems I there’s a section called “Solving The Nearest-Point-On-Curve Problem” where some public domain C code is given that computes the closest point on a cubic Bezier curve to a given input point. The C code is really old and actually calls malloc in a recursive function, presumably to avoid stack overflows on old machines with minimal memory. This can be converted to use a stack data structure on the process stack space, all without much work for the user.

After a few minutes I was able to set up a cool demo:


This tool would be great for some kind of interactive spline editor. It would also be perfect for colliding circles (or spheres in 3D), making the curve itself a potential candidate for a physics engine collider type. Since the curve is defined by control points, the curve can actually deform at run-time without any big complications.

Approximating the Curve

In reality a spline like this would probably just be “rasterized” to some kind of approximate shape that is easier for a physics engine to collide, rather than dealing with the exact spline itself. This is likely due to there not being a straightforward way to compute the closest point on the spline to a plane (at least, I sure don’t know how to write that code).

One tactic would be to interpolate a few points out of the curve (and by few I mean as many as needed), and then use these points as an approximate representation of the curve. This could be like converting the curve to a line strip, which might make interactivity much simpler in many cases.

Generating these points can be a piece of cake. Personally I’d just write out some lerps and call it good. If we look at the simplest Bezier curve with 2 control points it’s easy to see it’s just a single lerp (these animations were taken from the Wikipedia page for Bezier curves):



If we add a third point we can perform two lerps, and then lerp the results together. Intuitively this gives us a quadratic Bezier curve, and looks like:



Following suit we can use 4 control points and create a cubic Bezier curve:



An easy way to convert a cubic Bezier curve to a line strip involves sampling points on curve at time t, where t is on the interval of 0 to 1. Here’s some example code:

For example, Box2D defines a “chain shape” as a line strip of segments. Generally this strip is used to approximate curved, non-moving surfaces. Splines can very easily be converted to strips and sent to Box2D to define a chain shape.

Some Spliney Links


Minimum Bounding Sphere

Today I was reading this article by Ericson. The article is on the derivation involved in writing a function to calculate the minimum bounding sphere of a set of three or four points. The topic itself isn’t very interesting (to me) when compared to the algebra itself. Ericson skips the algebra I had trouble with, so my article here will serve the purpose to record the algebra he skipped (in case other readers had trouble, or I ever want to reference the math again). The article I am writing here directly follows Ericson’s article, so I will assume readers have read and understood Ericson’s article.

Say we have a set of three points S = {A, B, C}, which can be thought of as a triangle. The problem of calculating the minimum bounding sphere of S involves checking to see if the circumcenter of S lies within or outside the triangle S. Lets name the circumcenter P. One immediate implementation involves computing P, followed by the computation of the barycentric coordinates of P with respect to S. These coordinates can be used check if P lay within S.

However computing P of S followed by the barycentric coordinates of P with respect to S involves a lot of redundant computations. Furthermore, if P does not lay within S then P itself does not need to be computed at all, since only the barycentric coordinates were needed.

It would be nice to compute barycentric coordinates of P with respect to S directly, and only if necessary construct P.

P in relation to S can be defined as:

P = A + s*(B – A) + t*(C – A)

Where \(s\) and \(t\) are barycentric coordinates of S such that \(s – t – (1.0 – s – t) = 0\) if P is within S. Since P is equidistant from A, B and C, we can express P in relation to S with the following:

dot(P – B, P – B) = dot(P – A, P – A)

dot(P – C, P – C) = dot(P – A, P – A)

The interesting part involves plugging \eqref{eq1} into \eqref{eq2} and \eqref{eq3} followed by a collection of terms. Since I myself am a noob when it comes to algebra I had to go look up algebraic properties of the dot product to make sure I did not screw anything up. In particular these few rules are useful:

dot(u, v + w) = dot(u, v) + dot(u, w) \\
dot(s * u, v) = s*dot(u, v) \\
-dot(u – v, u – w) = dot(v – u, u – w)

Lets go over substituting \eqref{eq1} into \eqref{eq2} directly, however lets only consider the first part of \eqref{eq2} \(Dot(P – B, P – B)\) to keep it simple:

dot(A + s*(B – A) + t*(C – A) – B, A + s*(B – A) + t*(C – A) – B)

Since things immediately get really long some substitutions help to prevent errors:

u = A – B\\
v = s*(B – A) + t*(C – A)\\
w = u + v \\

Dot(w, u + v)

Dot(w, u) + Dot(w, v)

Dot(u + v, u) + Dot(u + v, v)

Dot(u, u) + Dot(v, u) + Dot(u, v) + Dot(v, v)

If I were to expand \eqref{eq10} on this webpage it would not fit. Instead we will make use of a few more substitutions and then arrive in the final form:

x = s*(B – A) \\
y = t*(C – A) \\
x + y = v

Dot(A – B, A – B) + 2*Dot(A – B, x + y) + Dot(x + y, x + y)

By following the same process we can finish the substitution and notice that:

Dot(P – A, P – A) = Dot(x + y, x + y)

The final form of the substitution would be:

Dot(A – B, A – B) + \\ 2*Dot(A – B, x + y) + \\ Dot(x + y, x + y) = Dot(x + y, x + y)

Dot(A – B, A – B) + 2*Dot(A – B, x + y) = 0

Dot(A – B, A – B) + \\ 2*Dot(A – B, s*(B – A)) + \\ 2*Dot(A – B, t*(C – A)) = 0

s*Dot(B – A, B – A) + t*Dot(B – A, C – A) = \\ (1/2)*Dot(B – A, B – A)

This final equation \eqref{eq17} matches exactly what Ericson came up with on his own blog. Through a similar process \eqref{eq1} can be substituted into \eqref{eq3}, which would result in:

s*Dot(C – A, B – A) + t*dot(C – A, C – A) = \\ (1/2)*dot(C – A, C – A)

Computing AABB Trick (Loop Trick)

Lately I noticed a small trick that applies to loops when trying to find a minimum or maximum of some values. Usually I just apply the trick to loops where I need to compute an AABB over some geometry. I think I noticed this trick when reading a for loop Erin Catto wrote in some of Box2D’s internal code.

The trick is super simple: just process the first element outside of the loop to set up your initial conditions, then form your loop to skip the first element. An assumption would be made that there’s at least one element in the array to process. Here’s an example for computing an AABB:

Usually I myself would have written this kind of code like so and not given any more thought to it:

This second code chunk is arguably just slightly more esoteric and is definitely a little less efficient for no good reason.

One could also skip the first element when finding the min/max of any sort of array, like for example: dot product results. Though simple it’s pretty nice to find small ways to write slightly better code.

Robust Parallel Vector Test

In an effort to try to find some robust tests for parallel vectors I’ve decided to create some slides talking about the solutions I’ve come across. Hopefully one day some reader can comment or email about some alternatives or better solutions! I decided to write this blog post in slide format since Microsoft’s Powerpoint is pretty snazzy for making quick diagrams and whatnot :)

Here’s the link to Ericson’s blog post on tolerances: Tolerances Revisited. Please note there’s a little bit of ambiguity about whether the test should care if the vectors point in the same direction or not. In general it doesn’t really matter since the point of the slides is numeric robustness.

The gist of the slides is that the scale of the vectors in question matters for certain applications. Computing a good relative epsilon seems difficult, and maybe different epsilon calculations would be good for different applications. I’m not sure!

Here’s a demo you can try out yourself to test two of the solutions from the slides (tests returning true should be considered false positives):

Output of the program is:

Download (PDF, Unknown)

Distance Point to Line Segment

2014-07-23 00.16.37

Understanding the dot product will enable one to come up with an algorithm to compute the distance between a point and line in two or three dimensions without much fuss. This is a math problem but discussion and focus is on writing a function to compute distance.

Dot Product Intuition

The dot product comes from the law of cosines. Here’s the formula:

c^2 = a^2 + b^2 – 2ab * cos\;\gamma

This is just an equation that relates the cosine of an angle within a triangle to its various side lengths a, b and c. The Wikipedia page (link above) does a nice job of explaining this. Equation \eqref{eq1} can be rewritten as:

c^2 – a^2 – b^2 = -2ab * cos\;\gamma

The right hand side equation \eqref{eq2} is interesting! Lets say that instead of writing the equation with side lengths a, b and c, it is written with two vectors: u and v. The third side can be represented as u - v. Re-writing equation \eqref{eq2} in vector notation yields:

|u\;-\;v|^2 – |u|^2 – |v|^2 = -2|u||v| * cos\;\gamma

Which can be expressed in scalar form as:

(u_x\;-\;v_x)^2 + (u_y\;-\;v_y)^2 + (u_z\;-\;v_z)^2\;- \\
(u_{x}^2\;+\;u_{y}^2\;+\;u_{z}^2) – (v_{x}^2\;+\;v_{y}^2\;+\;v_{z}^2)\;= \\
-2|u||v| * cos\;\gamma

By crossing out some redundant terms, and getting rid of the -2 on each side of the equation, this ugly equation can be turned into a much more approachable version:

u_x v_x + u_y v_y + u_w v_w = |u||v| * cos\;\gamma

Equation \eqref{eq5} is the equation for the dot product. If both u and v are unit vectors then the equation will simplify to:

dot(\;\hat{u},\;\hat{v}\;) = cos\;\gamma

If u and v are not unit vectors equation \eqref{eq5} says that the dot product between both vectors is equal to cos( γ ) that has been scaled by the lengths of u and v. This is a nice thing to know! For example: the squared length of a vector is just itself dotted with itself.

If u is a unit vector and v is not, then dot( u, v ) will return the distance in which v travels in the u direction. This is useful for understanding the plane equation in three dimensions (or any other dimension):


The normal of a plane would be the vector: { a, b, c }. If this normal is a unit vector, then d represents the distance to the plane from the origin. If the normal is not a unit vector then d is scaled by the length of the normal.

To compute the distance of a point to this plane any point can be substituted into the plane equation, assuming the normal of the plane equation is of unit length. This operation is computing the distance along the normal a given point travels. The subtraction by d can be viewed as “translating the plane to the origin” in order to convert the distance along the normal, to a distance to the plane.

Writing the Function: Distance Point to Line

The simplest function for computing the distance to a plane (or line in 2D) would be to place a point into the plane equation. This means that we’ll have to either compute the plane equation in 2D if all we have are two points to represent the plane, and in 3D find a new tactic altogether since planes in 3D are not lines.

In my own experience I’ve found it most common to have a line in the form of two points in order to represent the parametric equation of a line. Two points can come from a triangle, a mesh edge, or two pieces of world geometry.

To setup the problem lets outline the function to be created as so:

The two parameters a and b are used to define the line segment itself. The direction of the line would be the vector b – a.

After a brief visit to the Wikipedia page for this exact problem I quickly wrote down my own derivation of the formula they have on their page. Take a look at this image I drew:

2014-07-23 00.16.37

The problem of finding the distance of a point to a line makes use of finding the vector that points from p to the closest point on the line ab. From the above picture: a simple way to calculate this vector would be to subtract away the portion of a – p that travels along the vector ab.

The part of a – p that travels along ab can be calculated by projecting a – p onto ab. This projection is described in the previous section about the dot product intuition.

Given the vector d the distance from p to ab is just sqrt( dot( d, d ) ). The sqrt operation can be omit entirely to compute a distance squared. Our function may now look like:

This function is quite nice because it will never return a negative number. There is a popular version of this function that performs a division operation. Given a very small line segment as input for ab it is entirely possible to have the following function return a negative number:

It’s very misleading to have a function called “square distance” or “distance” to return a negative number. Passing in the result of this function to a sqrt function call can result in NaNs and be really nasty to deal with.

Barycentric Coordinates – Segments

A full discussion of barycentric coordinates is way out of scope here. However, they can be used to compute distance from a point to line segment. The segment portion of the code just clamps a point projected into the line within the bounds of a and b.

Assuming readers are a little more comfortable with the dot product than I was when I first started programming, the following function should make sense:

This function can be adapted pretty easily to compute the closest point on the line segment to p instead of returning a scalar. The idea is to use the vector from p to the closest position on ab to project p onto the segment ab.

The above function works by computing barycentric coordinates of p relative to ab. The coordinates are scaled by the length of ab so the second if statement must be adapted slightly. If the direction ab were normalized then the second if statement would be a comparison with the value 1, which should make sense for barycentric coordinates.

Sample Program

Here’s a sample program you can try out yourself:

The output is: “Distance squared: 2.117647″.

Deriving OBB to OBB Intersection and Manifold Generation


The Oriented Bounding Box (OBB) is commonly used all over the place in terms of simulation and game logic. Although often times OBBs aren’t actually bounding anything, and really just represent an implicitly defined box.

OBBs commonly come in 2D and 3D form, though there aren’t all that many resources on the internet describing how to derive the 3D case. The ability to write robust collision detection, and more importantly manifold generation (contact data needed to resolve the collision), relies on a good understanding of linear algebra and geometry.

Readers will need to know about affine and linear transformations, and have intuition about the dot and cross products. To learn these preliminaries I highly recommend the book “Essential Mathematics” by Jim Van Verth. I will also be assuming readers are familiar with the Separating Axis Theorem (SAT). If you aren’t familiar with SAT please do a short preliminary research on this topic.

The fastest test I’m aware of for detection collision between two OBBs, in either 2 or 3 dimensions, makes use of the Separating Axis Theorem. In both cases it is beneficial to transform one OBB into the space of the other in order to simplify particular calculations.

Usually an OBB can be defined something like this in C++:

6 Face Axes

As with collision detection against two AABBs, OBBs have x y and z axes that correspond to three different potential axes of separation. Since an OBB has local oriented axes they will most of the time not line up with the standard Euclidean basis. To simplify the problem of projecting OBB \(A\) onto the the axes of OBB \(B\), all computations ought to be done within the space of \(A\). This is just a convention, and the entire derivation could have been done in the space of \(B\), or even in world space (world space would have more complex, and slow, computations).

The translation vector \(t\) is defined as \(t = R_{a}^{T} * (C_b – C_a)\), where \(C_i\) denotes the center of an OBB, and \(R\) denotes the rotation of an OBB. \(t\) points from \(A\)’s center to \(B\)’s center, within the space of \(A\). The relative rotation matrix to transform from B’s frame to A’s frame is defined as \(C = R_{a}^{T} * R_b\). \(C\) can be used to transform \(B\)’s local axes into the frame of \(A\). In \(A\)’s own frame it’s local axes line up with the standard Euclidean basis, and can be directly computed with the half width extents \(e_a\).

To calculate the separation \(s\) along the local axis \(l\) from an OBB examine the following diagram and matching equation:


s = |t \cdot l| – (|a \cdot l| + |(C * b) \cdot l|)

Due to symmetry the absolute value of a projection can be used in order to gather a value that can be used to represent a projection along a given direction. In our case we aren’t really interested in the sign of the projection and always want the projection towards the other OBB. I understand that newer programmers may have some trouble translating math into code, so here’s an example of eq.\eqref{eq1} as C++:

s is the separation along an axis. t is the translation vector from center of one box to another. l is the axis we are projecting onto (in this case A’s x axis).

9 Edge Axes

In 2D eq.\eqref{eq1} is all that is needed to detect overlap upon any given possible axis of separation. In 3D a new topological feature arises on the surface of geometry: edges. Because edges form entirely separate Voronoi regions from face surfaces they can also provide possible axes of separation. With OBBs it is fast to test all possible unique edge axes of separation.

Given two edges upon two OBBs in 3D space the vector perpendicular to both would be the possible axis of separation, and can be computed with a cross product. A possible 9 axes of separation arise from crossing each local axis of \(A\) with the each local axis of \(B\).

Since \(C\) is composed of a concatenation between \(R_{a}^{T}\) and \(R_b\), the rows of this matrix can be thought of as the local axes of \(B\) in the space of \(A\). When a cross product between two axes is required to test separation, the axis can be pulled directly out of \(C\). This is much faster than computing a cross product in world space between two arbitrary axis vectors.

Lets derive a cross product between \(A\)’s x axis and \(B\)’s z axis, all within the space of \(A\). In \(A\)’s space \(A\)’s x axis is \(\begin{Bmatrix}1 & 0 & 0\end{Bmatrix}\). In \(A\)’s space \(B\)’s z axis is \(\begin{Bmatrix}C_{20} & C_{21} & C_{22}\end{Bmatrix}\). The axis \(l\) is just the cross product of these two axes. Since \(A\)’s x axis is composed largely of the number \(0\) elements of \(C\) can be used to directly create the result of the cross product, without actually running a cross product function. We arrive to the conclusion that for our current case \(l = \begin{Bmatrix}0 & -C_{22} & C_{21}\end{Bmatrix}\). This particular \(l\) axis can be used in eq.\eqref{eq1} to solve for overlap.

Since the matrix \(C\) is used to derive cross product results the sign of the elements will be dependent on the orientation of the two OBBs. Additionally, the cross product will negate on of the terms in the final axis. Since eq.\eqref{eq1} requires only absolute values it may be a good idea to precompute abs(\(C\)), which takes the absolute value element by element. This computation can be interleaved between the different axis checks for the first 6 checks, only computed as needed. This may let an early out prevent needless computation.

The same idea can be used to derive this axis in the space of \(B\) by crossing \(B\)’s local z axis with \(A\)’s local x axis transformed to the space of \(B\). To do this it may be easy to use \(C^T\), which represents \(R_{b}^{T} * A\). More information about this can be found in my previous article about the properties of transposing a matrix concatenation.

Eventually all axes can be derived on paper and translated into code.

Manifold Generation

A collision manifold would be a small bit of information about the nature of a collision. This information is usually used to resolve the collision during physical simulation. Once an axis of least penetration is found the manifold must be generated.

For SAT manifold generation is quite straight forward. The entirety of the details are expressed very well by Dirk Gregorius’ great online tutorial from GDC 2013. The slides are openly available for free, and currently hosted here generously by Erin Catto. Readers are encouraged to familiarize themselves with the idea of reference and incident faces, as well as the idea of reference face side planes. For a 2D analogue I have provided my own slides on this topic.

There are only two different cases of collision that should be handled: edge to edge and something to face. Vertex to edge, vertex to face, and vertex to vertex resolutions are ill-defined. On top of this numerical imprecision makes detecting such cases difficult. It makes more sense in every way to skip these three cases and stick to the original two.

Something to Face

Something has hit the face of one of the OBBs. It can be anywhere from 1 to 4 vertices from the face of the opposing OBB. The entire intersecting volume between two OBBs does not need to be considered due to the nature of SAT and the axis of minimum penetration. Only two faces from each OBB need be considered for this case, and are by convention named the reference and incident faces, one from each OBB.

Finding the incident face given an axis of separation can be found in constant time with implicit geometry like an OBB; the problem reduces to that of finding a vector that most faces a given vector, or in other words the most negative dot product.

Once the incident face is identified it is clipped to the side planes of the reference face. There may be some tricks involving basis transformations here to simplify the clipping routine. In order to perform clipping the incident face vertices must be calculated. Given a face normal from an OBB all four face vertices can be computed with the normal and the half width extents of the OBB. In similar fashion the side planes of the reference face can also be directly computed.

Edge to Edge

The edge to edge case requires a clever constant-time computation of a supporting edge. A supporting edge is the edge most pointed to by a given axis of separation. The computation of this edge may be tricky to derive yourself, though there do exist online references. I suggest checking out Open Dynamics Engine (ODE) by Russell Smith, or the Cyclone Engine by Ian Millington for information on how to compute a supporting edge.

Once both supporting edges are calculated the closest two points between each edge can be calculated and used as a single contact point of which to resolve upon.

Numeric Stability

In all uses of OBB to OBB intersection cross products with nearly parallel vectors pose a problem. The best solution is to detect when any two axes from each OBB are parallel and skip all nine cross product axes. It can be visually shown that no false positives will come from skipping the cross product axes in this case.

A simpler solution, as proposed by Gottschalk ’00 in his paper “Collision Queries using Oriented Bounding Boxes”, which made its way into Ericson’s book “Real-Time Collision Detection”, is to bias all elements of \(C\) with a small epsilon to drive near-zero axes away from a degenerate case. For other axes the epsilon value should be tuned such that it doesn’t have much an impact upon results.

If using this collision routine to resolve collision during physical simulation certain tuning my be appropriate. If an iterative resolver is used it may be preferred to have slightly more consistent manifold generation, as opposed to exact manifold generation. For details I suggest researching works by Erin Catto and comments made by others around the Bullet Forums.


Hopefully by now readers realize that this whole time we’ve been simplifying the problem of OBB to OBB to the problem of AABB to OBB. When an actual AABB needs to be collided against an OBB simplifications can occur. Many simplifications revolve around the non-necessity to transform things into the basis of the AABB, since world coordinates should suffice.

One Axis Aligned (or Two?)

Some games have OBBs that only orient on some of their axes instead of all of them. One can immediately realize that no cross product checks need to be performed in any case, if two OBBs have orientation locked on the same axis. Other optimizations will occur as well, making some axis tests reduced to that of AABB overlap testing.


Collision detection is difficult. Writing robust collision detection code requires a good mathematical foundation as well as the ability to write efficient code. This slight mixture of fields of study may be difficult to learn all at once, but it all gets easier with time and practice and the rewards can be very fulfilling.

I’d like to thank my friend Nathan Carlson for teaching much of this information to me, and for his insightful discussions.

Here is a finished implementation: link to github file.

Buoyancy, Rigid Bodies and Water Surface


I’ve spent the last couple weeks researching eigenvalue decomposition and solving cubic polynomials in order to simulate a liquid surface and polyhedral buoyancy! I gave a lecture about this topic at my university and hope the slides can be of interest or help to someone. I’ve attached a couple demo GIFs to the bottom of this post.

Download (PDF, Unknown)



Dynamically Slicing Shapes

Just this morning I finished up a small open source project, along with a huge article at gamedev.tutsplus. The open source project is an implementation of Sutherland-Hodgman clipping which can be used in many advanced shape and mesh manipulation techniques.


Gif of the open source project in action. Slicing a quad interactively.