Convex hull algorithms

In this post we’re going to go over ways we can build convex shapes from a set of arbitrary points. This process is done by finding the points that connect to encapsulate all points in a given set (whilst remaining convex), we call this a convex hull. In game development convex hulls are commonly used as optimizations for collision. Usually perfect collision detection for high triangle meshes would be resource hogs. Instead we tend to calculate their convex hulls so that we can do fast and somewhat accurate collision tests on them. We don’t use the convex hull for every mesh, but most static objects like chairs or tables, etc. probably are using convex hulls for collision.

This post will go over two algorithms both in 2D. The first is limited specifically to 2D, but it will help us understand the basic idea behind what a convex hull algorithm is achieving. The second is more robust as it can be expanded to 3D although that will not be covered in the post as it is significantly more complex (perhaps that is a post for another time).

---

Prerequisites

• Dot product
• Cross product
• Vector math

Andrew’s monotone convex hull

This algorithm works by sorting the points within the set from left to right and then building the lower hull, then the upper hull, and putting them together. The upper hull (blue) simply refers to the top half of the convex hull and the lower hull (red) refers to the bottom half of the polygon.

Sorting the points from left to right is very easy. We can just use the table.sort function alongside vector projections onto the x-axis to put the points in order of left to right.

function dot(a, b)
    return a.x * b.x + a.y * b.y;
end;

function monotoneHull(set)
    -- sort points from left to right
    table.sort(set, function(a, b) return dot(a, Vector2.new(1, 0)) < dot(b, Vector2.new(1, 0)); end);
end;

Now that we have the points in order from left to right we’ll start with the lower hull. In order to find which points should be a part of the hull and which shouldn’t we’re going to need to implement something we’ll call the 2D cross product. Normally the cross product is a thing we can only use in 3D, but for our purposes we’ll define it as the z component alone.

function cross(a, b)
    -- z-component of normal cross product
    return a.x * b.y - a.y * b.x;
end;

The reason this value is important to us is because if we were to pretend all our vector2 values were vector3 values with a z-component of zero, then when we crossed them we each other we’d either get a vector pointing out of our screen or pointing into our screen. This means that we’d end up getting a vector 3 value that had zeros as both its x and y components and either a positive or negative z-value. If we understand the right-hand rule then whether that z component is positive or negative tells us a lot about the two vectors in question, specifically if one is above or below the other relative to a third vector.

To further show what’s meant by this let’s use the above plotted points to go over a few iterations of the monotone hull algorithm trying to find the lower hull. To start we’ll pick the first two points from left to right and subtract from the second from the first. We’ll then take the third point and also subtract it from the first. We’ll assume these are all part of the lower hull for the time being. This will leave us with two vectors that are relative to the same origin (u and v).

We can then take the cross product of u against v which in this case will give us a positive value (pointing out of the screen if we did the 3D cross product). Since we got a positive value we know that the third point added is above the second point. Since we’re building the lower hull you might think we should remove the third point and continue testing forward. The problem with doing this is that just because a point is below is the previous one doesn’t mean it’s not part of the lower hull, imagine if our set was only made of these three points! As such we need to do further testing.

Let’s find the next point from left to right (the fourth point) and take that point and the third point and subtract them from the second point giving us new u and v vectors.

Once again we can cross u onto v which gives us a positive value. This means we repeat the same process as we did for the last set of points. We find the next point (in this case the fifth) and we subtract it from the third point to get our new v vector. We then take the fourth point and subtract it from the third point to get our new u vector.

Aha! Now when we cross u against v we get a negative number! This means that without a doubt there is a point to the right side of point four that is lower! This means it’s not possible for point four to be part of the lower hull and as such we can remove it without a doubt. This leaves us with a hull of point 1, point 2, point 3, and point 5 so far. Now that we’ve removed a point we need to go back to point 3 and check if point 5 is also lower which would mean point 3 is also not part of the lower hull. We do this by once again subtracting point 3 from point 2 for our u vector and subtracting point 5 from point 2 for our v vector. Sure enough when we cross u against v we get a negative value meaning we can remove point 3 from the hull.

Once again we must repeat the process and check if point 2 is lower than point 5. To do this we subtract point 2 from point 1 for the u vector and we subtract point 5 from point 1 for v vector. When we cross these vectors we get a positive value meaning that point 5 is higher than point 2 which means point 2 is also part of the lower hull.

Since we’ve gone back and checked all the points we previously had in the hull we now move forward by getting the next point in the set (point 6) and doing the check as we were doing above. Once we’ve added and checked all the points in the set we’re left with the points in order from left to right on the lower convex hull!

function monotoneHull(set)
    table.sort(set, function(a, b) return dot(a, Vector2.new(1, 0)) < dot(b, Vector2.new(1, 0)); end);
    -- lower hull
    local low = {}; -- start with empty lower hull
    for i = 1, #set do
        local p = set[i]; -- get next point
        -- make sure we have at least two points in addition to the next point
        while #low >= 2 and cross(low[#low] - low[#low-1], p - low[#low-1]) <= 0 do
            -- if 2d cross was equal to zero or negative remove the last point currently in the lower hull
            table.remove(low, #low);
        end;
        -- add point to lower hull
        table.insert(low, p);
    end;
    return low;
end;

To make sure that everything’s working let’s plug in the points used in the graph above and visually check that the points returned are in order and make up the lower hull.

local test = {
    Vector2.new(1204, 734);
    Vector2.new(969, 385);
    Vector2.new(599, 304);
    Vector2.new(701, 228);
    Vector2.new(726, 428);
    Vector2.new(1347, 334);
    Vector2.new(164, 5);
    Vector2.new(828, 126);
    Vector2.new(14, 284);
    Vector2.new(726, 428);
    Vector2.new(701, 228);
    Vector2.new(599, 304);
    Vector2.new(498, 112);
    Vector2.new(228, 738);
    Vector2.new(202, 126);
};

for _, p in next, monotoneHull(test) do
    print(p);
end;

-- Output:
> 14, 284
> 164, 5
> 828, 126
> 1347, 334

Awesome! Everything checks out for the lower hull!

The next question is how do we get the points in the upper hull? Luckily for us the answer is pretty much the exact same as it is for the lower hull. The only difference is that this time we travel across our set from right to left, other than that we don’t have to change anything because in traveling the opposite direction our cross values will also swap, meaning we should still remove a point if we get a negative cross or keep it if we get a positive.

The one thing we have to be aware of is that both the upper and lower hull with both have the furthest most left and right points. As such when we take the two parts of the full hull we have to remember to remove the left and right from at least one of the upper or lower hull else our final hull will have duplicate points.

function monotoneHull(set)
    table.sort(set, function(a, b) return dot(a, Vector2.new(1, 0)) < dot(b, Vec-tor2.new(1, 0)); end);
    -- lower hull
    local low = {};
    for i = 1, #set do
        local p = set[i];
        while #low >= 2 and cross(low[#low] - low[#low-1], p - low[#low-1]) <= 0 do
            table.remove(low, #low);
        end;
        table.insert(low, p);
    end;
    -- upper hull
    local up = {};
    for i = #set, 1, -1 do
        local p = set[i];
        while #up >= 2 and cross(up[#up] - up[#up-1], p - up[#up-1]) <= 0 do
            table.remove(up, #up);
        end;
        table.insert(up, p);
    end;
    -- account for duplicates
    table.remove(low, #low);
    table.remove(up, #up);
    -- put the hulls together
    for _, p in next, up do table.insert(low, p); end;
    -- return the hull
    return low;
end;

Sure enough when we run this code we get the convex hull!

Quick hull 2D

The quick hull algorithm is all about dividing and conquering. Throughout the process our goal is going to be to build a simplex (in 2D this means a triangle) and use the edge normals in addition to the process of elimination to get the points that make up the convex hull. Just like last time we’ll go over a few iterations so that we can fully understand what quick hull is doing.

To start off we want to pick an arbitrary direction that we can sort the points based on. In the previous algorithm we used (1, 0) to sort from left to right. With quick hull we can pick any initial direction we want, however for simplicity sake we’re just going to use (1, 0) again as it’s pretty straightforward to wrap your head around when it comes to scalar projections.

Once we’ve picked that initial direction we want to get the two most extreme points projected along it. We already learned how to do with the dot product and the table.sort function in the previous algorithm, so this will be a cinch!

function quickHull(set)
    local d = Vector2.new(1, 0);
    table.sort(set, function(a, b) return dot(a, d) > dot(b, d); end);
    local p0, p1 = set[1], set[#set];
end;

Because both those points are the extremes in a given direction that we’re looking at for the first time we know they’re both a part of the convex hull. Our next step is going to be to find the vector between them and set the new direction as perpendicular to that vector. Since there are two possible perpendicular vectors you might be wondering which one to pick? Luckily for us it doesn’t matter because we’re going to come back and use the unused direction later. For example sake we’ll pick the perpendicular vector going “downwards”.

We can then use this new d vector to get the point that is furthest in that direction thereby allowing us to create a simplex. We can then take this new point and order so that it represents a connection between the first and second points.

Once we’re in a position where we have a triangle the process becomes very mechanical. We go from left to right taking the perpendicular unit vectors facing outwards to each newly created edge and find any extreme points in the set that aren’t already part of the hull and lie beyond the edge.

In the case of the v vector when we find the point in the set that is furthest in that direction we get a point that’s already in our hull. As such move onto the next edge vector, u. When we test u we find that the furthest point in that direction that isn’t part of the hull and lies beyond the edge we’re testing against we get (828, 126). As such we add it to our convex hull as a connection between (164, 5) and (1347, 334).

We then once again go from left to right and check if there’s an extreme point in the direction of v past the edge’s outward normal it represents. In this case no points in the set pass that test and as a result we move onto vector u. When we do the test on vector u we also find there isn’t an extreme point in that direction past the edge either.

So we’re now at the point where you might be wondering what we test next? We checked all the edges in the simplex and built one half of the hull, but what about the points higher up? This is where we now take the other perpendicular vector to the initial two points in our hull. Once we’ve set that in motion we just repeat the exact same process as above by checking the simplex normals the only difference is that we now check the extremes from right to left.

Before we can get to coding this there’s two things we’ve neglected to talk about and that’s how we get the edge normals once we have a simplex, and how to check if an extreme is past an edge or not.

Vector triple product

When it comes to the edge normal, we can’t just randomly pick one of two possible perpendicular vectors we specifically need the one facing outwards. As such we need a way to approach this problem accurately and without guesswork. Let’s look at a simple test case with only a few points so that we don’t over complicate things.

If we take u and cross it with v we will get a vector pointing out of the screen. If we then take that vector and cross it with u again we get a vector that’s perpendicular to u and is facing towards C. If we take (v x u) x u then we get the perpendicular vector facing away from C.

With this in mind we can apply it to our simplex to get the proper edge facing normal. Let’s look at the first simplex we had in the example above. We have three points, if we take the second and subtract it from the first we get vector v, likewise if we take the third and subtract it from the first we get vector u.

If we take we take (v x u) x u and normalize it we get the vector perpendicular to v facing outwards! Awesome! We have a way to solve for edge normals.

function tripleProduct(a, b, c)
    -- (a x b) x c
    -- convert to vector3 temporarily
    local a3 = Vector3.new(a.x, a.y, 0);
    local b3 = Vector3.new(b.x, b.y, 0);
    local c3 = Vector3.new(c.x, c.y, 0);
    local trip = a3:Cross(b3):Cross(c3);
    -- convert back to vector2
    return Vector2.new(trip.x, trip.y);
end;

The above code will work fine, but you might not like the idea of having to convert your vector2 values into 3D, doing the math, and then converting back. Luckily for us there’s an optional expansion that allows us to do the same operation with only some simple vector multiplication, subtraction, and the dot product thereby allowing us to stay in the realm of 2D.

Optional expansion

To start let’s expand (A x B) x C and see where than leaves us:

Next to help us expand we’re going to add and subtract the product of all three vectors to each component. This doesn’t actually change anything because we’re subtracting exactly what we’re adding, but it helps us simplify the above into the expanded form.

Now we can take each component and factor some stuff out

And would you look at that, we’re left with the dot product and some simple vector multiplication! We can write the above form as:

function tripleProduct(a, b, c)
    -- (a x b) x c = b * (c . a) - a * (c . b);
    return b * dot(c, a) - a * dot(c, b);
end;

Checking if a point is beyond an edge

The second thing we need is a way to check if a point is beyond an edge or not. Luckily this process is very easy with one of our simple dot product definitions:

If we take the extreme point relative to one of the edge vertices and dot it with the edge normal we know that since magnitudes are always positive if the theta is greater than 90 degree we’ll get a negative value (or zero) and if it’s less than or equal to 90 degrees we’ll get a positive value. With that in mind if we can use this property to tell if a point lies beyond an edge or not.

Putting quickhull into code

Alright, we now have all the pieces we need to make the quick hull algorithm work. So taking what we’ve learned we can put the above into code

function quickHull(set)
    -- initial direction to get two points
    local d = Vector2.new(1, 0);
    table.sort(set, function(a, b) return dot(a, d) > dot(b, d); end);
    local p0, p1 = set[1], set[#set]; -- the max left and right points
     -- two seperate tables one for the hull and one for points that are in the hull be we need to help us calculate things
    local points, hull = {p0, p1}, {};
    -- get the vector that splits the set into two
    local originv, firstpass = p1 - p0, true;
    -- get one of the perpendicular directions.
    d = Vector2.new(originv.y, -originv.x);
    while #points > 1 do
        -- get the extreme point in that direction
        table.sort(set, function(a, b) return dot(a, d) > dot(b, d); end);
        -- get the simplex points realtive to a single point for use in triple product later
        local furthest = set[1] - points[1];
        local edge = points[2] - points[1];
        -- see if the extreme point lies beyond the edge we're checking
        if furthest ~= edge and dot(furthest, d) > 0 then
            -- put the new point inbetween the first and second point
            points = {points[1], set[1], select(2, unpack(points))};
            -- find the edge normal
            d = tripleProduct(edge, furthest, furthest).unit;
        else -- no points lie beyond the edge
            local hullPoint = points[1]; -- we don't need this point anymore for calculations
            table.insert(hull, hullPoint); -- add it to the hull
            table.remove(points, 1); -- remove it from the points so we move to the next edge
            if #points > 1 then -- make sure there are more edges to check
                furthest = hullPoint - points[1];
                edge = points[2] - points[1];
                -- find the next edge direction
                d = tripleProduct(furthest, edge, edge).unit;
            elseif firstpass then -- if we removed all the points then we finished that side of the originv
                firstpass = false;
                points = {p1, p0}; -- now checking in reverse (left to right vs right to left)
                d = Vector2.new(-originv.y, originv.x).unit; -- update the direction
            end;
        end;
    end;
    -- return the calculated hull
    return hull;
end;

Conclusion

So that concludes two possible algorithms for making 2D convex hulls. Perhaps in the future we'll go over the broad strokes of 3D as any implementation will be hundreds of lines of code. As always I hope you enjoyed and learned something new. Thanks for reading!