|Previous - Triangles||Next - Boxes|
A polygon is a closed shape defined by three or more points. Polygons, like triangles, have a winding order defined by the order of the points. Unlike triangles however, polygons come in two additional flavors, convex or concave. Convex polygons, to which triangles and rectangles belong, have the property that if you take any two points belonging to the polygon, that all points on the line between those points also belong to the polygon. For a star polygon, which is convex, this is not the case for example.
Polygon area and winding order
The winding order of a polygon can be obtained by calculating a sum of cross products, which give us double the signed area of the polygon. The formula you find under the name shoelace formula or Gauss’s area formula usually comes under the following form
and gives us the area of the polygon. If we want the winding order, we don’t need the absolute value or division, thus we can look at
Which, if we switch some terms becomes
Which, after grouping some terms becomes
Which is is clearly several cross products
First of all, why the separate cross product? Well, we take two points in sequence every time, and at the end we need to wrap around, as the last point uses the first.
Next, how does this formula differ with the one we got for triangles? At first glance, it looks familiar but different. Both use cross products. Where the triangle formula used just one cross product, this one will use three for a triangle. And where the triangle formula used vectors originating from one of the points towards the remaining points, this formula uses vectors originating from the origin towards the points.
But do we get the same result for a triangle? Well let’s calculate.
Our triangle formula used to get the winding order was . So if A, B and C are written as , and we get that
Calculating the cross product gives us
Writing out the products gives us
The and cancel each other out so we get
Now let’s look at the new formula. For a 3 point polygon like our triangle we get
Which is exactly the same as our previous formula.
Now, the formula we are going to use in our code has yet another form
Where n+1 wraps back to 1. When we write out this formula for our triangle we get
Writing this out gives
Eliminating opposing pairs gives us
Which gives our original formula again. This form is used most of the time as it takes an addition, subtraction and multiplication rather than two multiplications and a subtraction, which is generally cheaper to do. Though if you use SIMD instructions, cross products is most likely the better choice.
function polygonWindingOrder(points) if #points < 6 then return nil end local sum = (points + points[#points-1]) * (points - points[#points]) for i=1, #points-3, 2 do sum = sum + (points[i+2] + points[i]) * (points[i+3] - points[i+1]) end if sum < 0 then return "ccw" else return "cw" end end function polygonArea(points) if #points < 6 then return nil end local sum = (points + points[#points-1]) * (points - points[#points]) for i=1, #points-3, 2 do sum = sum + (points[i+2] + points[i]) * (points[i+3] - points[i+1]) end return math.abs(sum) * 0.5; end
Convex or concave
Another important property of polygons is whether they are convex or concave. A convex polygon has the property that if you take any two points, the line between those points falls entirely within the polygon. Thus a rectangle or hexagon is convex, while a star polygon is not.
Whether a polygon is convex or concave has implications in use. Any convex polygon can for example be much easier split into triangles, while concave polygons which are not stars need more complex steps. The same goes for collision tests.
So how can we know whether a given polygon is convex or concave? Two adjacent sides of a convex polygon will always turn towards the polygon, they go around the polygon, not alternating towards and away like a star for example. Thus the angle between each two sides can tell us whether the polygon is convex or not. Instead of the angle we only need to know the sign of the sine however, as we only need to know whether all angles are greater or smaller than 0, which maps to a positive or negative sine. Since we are only interested in the sign and can disregard the magnitude, we can use the cross product.
So our final algorithm calculates the cross product between every two adjacent sides. We can exit early if the sign switches as we know the polygon will be concave in that case. If we pass the whole polygon with the same sign, we have a convex polygon.
function polygonType(points) if #points < 6 then return nil end local v1x, v1y = points[#points-1] - points, points[#points] - points local v2x, v2y = points - points, points - points local sign = v1x * v2y - v1y * v2x for i=1, #points-5, 2 do v1x, v1y = -v2x, -v2y v2x, v2y = points[i+4] - points[i+2], points[i+5] - points[i+3] if (sign * (v1x * v2y - v1y * v2x)) < 0 then return "concave" end end return "convex" end
Drag the vertices below to see how the cross product of adjacent edges changes depending on whether the vertex makes the polygon convex or concave. Green is used for a negative cross product, while red is used for a positive cross product.
Just like the triangle center, we take the sum of all points and divide by how many there are.
function polygonCenter(points) local centerX, centerY = points, points for i=3, #points, 2 do centerX = centerX + points[i] centerY = centerY + points[i+1] end local count = #points / 2 centerX = centerX / count centerY = centerY / count return centerX, centerY end
Point in polygon
One of the ways to check whether a point lies within a polygon is using the same method as we used with triangles, we look at the sign of cross products. However that are a lot of cross products when our polygon contains more than a few points. Of course, if we know the winding order and that the polygon is convex, we can reduce the amount of cross products. But since we want a general solution, we’re going to look at a faster way of doing it for all polygons. This method is courtesy of W. Randolph Franklin.
We still do a check for every edge, but we only look for edges which cross our y coordinate. When they do, we find the x coordinate on the edge and check whether it is larger than our x. The idea is that we count how many crossings of edges lie on the right side. If it is odd, we are inside the polygon, if it is even, we are outside the polygon.
The test to look whether y lies between y1 and y2, is written as (y1 >= y) != (y2 >= y). This test only passes if one of them is true, not both. So if (y1 >= y and y2 < y) or (y1 < y and y2 >= y). Moreover, note that it doesn’t pass when y1 == y2, which is important for the next test in order to avoid a division by zero.
To obtain x’ of the intersection at y’=y, to test whether x is smaller than x’, we can use the vector equation of the line.
Since we know y’ is equal to y, we can get t
Filling this t into the equation for x’ we get
function Polygon:contains(x, y) if #self.points < 6 then return false end local x1, y1, x2, y2 = self.points[#self.points-1], self.points[#self.points], self.points, self.points local c = false for i = 1, #self.points-1, 2 do if ((y1 >= y) ~= (y2 >= y)) and (x <= (x2 - x1) * (y - y1) / (y2 - y1) + x1) then c = not c end x1, y1 = x2, y2 x2, y2 = self.points[i+2], self.points[i+3] end return c end
You can see the result below. Blue is the polygon, red is the point that is being tested. Purple are points for which x is greater than the x of the red point. Green are points for which x is smaller than the x of the red point. When the amount of purple points is odd, the red point lies in the blue polygon. You can move both the red point and the polygon’s vertices to see what the result is for different setups.
Polygon collision detection
The following collision detection method only works with convex polygons. So if you have concave polygons, you either need to triangulate them or split them into convex polygons.
The collision detection algorithm we use, called SAT or separating axis theorem, is quite simple. Remember that we said that we could use vector projection to do collision detection? What we do is project both polygons onto a line. If the resulting line segments don’t overlap, the polygons don’t collide.
Of course if they do overlap, it doesn’t mean that they collide, there might be another line on which the projected polygons don’t overlap.
Luckily there’s a limit to the amount of projections we have to do. We only need to project onto the distinct normals of the polygons.
So the algorithm is basically
function projectPolygon(points, x, y, nx, ny) local min = (points-x)*nx+(points-y)*ny local max = min local s for i = 3, #points, 2 do s = (points[i]-x)*nx+(points[i+1]-y)*ny if s < min then min = s end if s > max then max = s end end return min, max end function polygonCollidesWith(points1, points2) -- We project both polygons on all normals of polygon 1 local x1, y1 = points1[#points1-1], points1[#points1] local x2, y2 = points1, points1 local nx, ny, min1, max1, min2, max2 -- For every edge of points1 for i = 1, #points1, 2 do -- Get edge normal by rotating the vector [x2-x1, y2-y1] by 90 degrees nx, ny = y2-y1, x1-x2; -- Project both polygons onto the normal min1, max1 = projectPolygon(points1, x2, y1, nx, ny) min2, max2 = projectPolygon(points2, x2, y1, nx, ny) -- If there is no overlap between the ranges, there is no collision if max1 <= min2 or min1 >= max2 then return false end x1, y1 = x2, y2 x2, y2 = points1[i+2], points1[i+3] end return true; end function polygonsCollide(points1, points2) return polygonCollidesWith(points1, points2) and polygonCollidesWith(points2, points1) end
Note that we don’t need a vector projection, a scalar projection is sufficient. The scalar projection of a point onto a vector gives us the position of the projected point as a one dimensional coordinate on the line defined by the vector. Note that we don’t even need to divide by the dot product of the vector we project onto as a scale will not change whether the ranges overlap or not.
So we take the unscaled scalar projection of each point in the polygon, and we only keep the minimum and maximum coordinate of these projections. This gives us our polygon flattened into a one dimensional range, or a line segment if we would project it back to 2D.
We do this for both polygons and check whether the ranges overlap, if they don’t, the polygons don’t collide. If they do, we check the next normal.
If we can’t find a normal for which the projections overlap, the polygons don’t collide.
You can test the algorithm below, by dragging vertices, edges or polygons. Note though that both your polygons should be convex.
While there is use for checking collisions without resolving the collision, think for example dragging and dropping a shape, like a puzzle piece, onto a hole, in many cases we do need to restore the polygons to a state where they don’t collide.
The idea is to use the overlap we found previously to displace one of the polygons. This means we need to keep the smallest overlap (since we want the smallest correction possible) as well as the direction vector in order to do the displacement.
Finally, we need to check whether our displacement vector points in the correct direction. Given that we might have two normals giving the same orientation but the opposite sign, a square for example, we need to check for this. We can use the dot product here. If the dot product of the displacement vector and the vector indicating the direction in which the displaced polygon lies with respect to the static polygon is smaller than 0, the vectors point in the opposite direction and we need to flip the displacement vector.
How do we determine the vector pointing from the static polygon to the displaced polygon? We take the centers of both polygons, and we subtract the center of the static polygon from the center of the displaced polygon.
function polygonCollidesWithV(points1, points2) -- We project both polygons on all normals of polygon 1 local x1, y1 = points1[#points1-1], points1[#points1]; local x2, y2 = points1, points1; local nx, ny, min1, max1, min2, max2; local overlap, smallestOverlap, length, overlapX, overlapY; -- For every edge of points1 for i = 1, #points1, 2 do -- Get edge normal by rotating the vector [x2-x1, y2-y1] by 90 degrees nx, ny = y2-y1, x1-x2 -- Project both polygons onto the normal min1, max1 = projectPolygon(points1, x2, y1, nx, ny); min2, max2 = projectPolygon(points2, x2, y1, nx, ny); -- If there is no overlap between the ranges, there is no collision if max1 <= min2 or min1 >= max2 then return false end overlap = math.min(max1, max2) - math.max(min1, min2); -- Our scalar projection wasn't scaled yet length = nx*nx+ny*ny; overlap = overlap/ length; -- And we need to take the length of the vector into account for the overlap length length = math.sqrt(length); -- Record smallest overlap if not smallestOverlap or overlap * length < smallestOverlap then smallestOverlap = overlap * length; overlapX = nx*overlap; overlapY = ny*overlap; end x1, y1 = x2, y2; x2, y2 = points1[i+2], points1[i+3]; end return smallestOverlap, overlapX, overlapY; end function polygonsCollideV(points1, points2) local overlap1, x1, y1 = polygonCollidesWithV(points1, points2) local overlap2, x2, y2 = polygonCollidesWithV(points2, points1) if overlap1 and overlap2 then if overlap1 < overlap2 then return overlap1, x1, y1 else return overlap2, x2, y2 end else return false end end
You can see how it works below. The green line is the vector from the center of the static polygon towards the displaced polygon. The red line is the vector which will move the displaced polygon away from the static polygon.
What if both polygons are moving? Then the displacement will be distributed to both polygons according to their mass. This will be explained once we get to physics.
While being able to detect collisions between convex polygons is nice, we often meet circles in games. Wheels, rolling stones or characters, they are all very round. Now, while a circle is most of the time drawn as a polygon approximating the actual circle, we don’t really need to use such an expensive polygon to test with. We can use approximately the same algorithm as before.
Projecting the circle is very simple, we only project the center point, and subtract and add the radius to get the minimum and maximum.
But what normal should we use? There are an infinite number of normals for a circle, and still far too many in an approximation.
We are in luck again, as the only axis we need to test for is the one defined by the center of the circle and the closest point in the polygon.
Note that unlike before, we normalize the normals, this is because otherwise we need to scale the radius of the circle into normal space, to be able to compare it with the other projections in normal space. To avoid this, we use a normalized vector, which doesn’t cost us more, as the alternative is to multiply the radius with the squared length of the vector and dividing it by the length of the vector.
function projectCircle(cx, cy, radius, x, y, nx, ny) local s = (cx-x)*nx + (cy-y)*ny; return s - radius, s + radius; end function polygonCircleCollideV(points, cx, cy, radius) -- We project both polygons on all normals of polygon 1 local x1, y1 = points[#points-1], points[#points] local x2, y2 = points, points local nx, ny, min1, max1, min2, max2 local overlap, smallestOverlap, length, overlapX, overlapY -- For every edge of points for i = 1, #points, 2 do -- Get edge normal by rotating the vector [x2-x1, y2-y1] by 90 degrees nx, ny = y2-y1, x1-x2 length = nx*nx+ny*ny length = math.sqrt(length); nx = nx / length ny = ny / length -- Project both polygons onto the normal min1, max1 = projectPolygon(points, x2, y1, nx, ny) min2, max2 = projectCircle(cx, cy, radius, x2, y1, nx, ny) -- If there is no overlap between the ranges, there is no collision if max1 <= min2 or min1 >= max2 then return false end overlap = math.min(max1, max2) - math.max(min1, min2) -- Record smallest overlap if not smallestOverlap or overlap < smallestOverlap then smallestOverlap = overlap overlapX = nx*overlap overlapY = ny*overlap end x1, y1 = x2, y2 x2, y2 = points[i+2], points[i+3] end -- For the axis from the circle center to the nearest point local index = 1 local nearest = cx*points+cy*points local distance for i = 3, #points, 2 do distance = cx*points[i]+cy*points[i+1] if distance < nearest then index = i nearest = distance end end local px, py = points[index], points[index+1] nx, ny = cx-px, cy-py length = nx*nx+ny*ny length = math.sqrt(length) nx = nx / length ny = ny / length -- Project both polygons onto the normal min1, max1 = projectPolygon(points, px, py, nx, ny) min2, max2 = projectCircle(cx, cy, radius, px, py, nx, ny) -- If there is no overlap between the ranges, there is no collision if max1 <= min2 or min1 >= max2 then return false end overlap = math.min(max1, max2) - math.max(min1, min2) -- Record smallest overlap if overlap < smallestOverlap then smallestOverlap = overlap overlapX = nx*overlap overlapY = ny*overlap end return smallestOverlap, overlapX, overlapY end
You can test the result below by dragging the circle towards the polygon. The additional circle axis comes into play when the circle is close to one of the vertices. Not that if the circle is too deep into the polygon, it needs several steps to get pushed out of it. This is because the maximum overlap, and thus displacement, will be the diameter of the circle.