Title: Determine whether a point is inside a polygon in Python

[This program draws points green if they're outside the polygon and red if they're inside it]

To use the program, left-click to select the polygon's points. Right-click to finish the polygon. Then move the mouse around the window and press the spacebar to make the program decide whether the mouse is inside or outside of the polygon. The program draws points green if they're outside the polygon and red if they're inside it.

One way to determine whether a point lies within a polygon is to add up the angles between the point and adjacent points on the polygon taken in order. For example, if the point in question is P and points A and B are adjacent on the polygon, then you look at the angle APB.

If the total of all the angles is zero, then the point is outside the polygon. If it's something else, the point is inside the polygon.

You can verify this intuitively if you think about a triangle or rectangle. If the point is outside of the shape, the clockwise angles are offset by counterclockwise angles.

If the point is inside the shape, then the angles continue sweeping in the same direction until you've gone full circle so the angles add up to ±360° or ±2π radians. If the shape is more complex, like the one in the picture at the top of this post, the angles may switch between clockwise and counterclockwise, but the result will be the same: ±2π.

If you look at the picture at the top of the post, you'll see that the point under the crosshair cursor is green and the sum of the angles is zero because the point is outside of the polygon.

Adding Angles

The following point_in_polygon function adds up the angles between a point and a polygon's vertices and returns True if the point is inside the polygon. It also returns the sum of the angles.

def point_in_polygon(points, pt): '''Return True if the point is in the polygon.''' # Add the angles from the point to each pair of vertices. total_angle = 0 num_points = len(points) for i in range(num_points): j = (i + 1) % num_points total_angle += get_radians(points[i], pt, points[j]) # If the total is near zero, the point is outside the polygon. return not math.isclose(total_angle, 0, abs_tol=0.000001), total_angle

This code initializes total_angle to 0. It then loops through the polygon's points, calls the get_radians function (described shortly) to get the angle between the point and two of the polygon's adjacent vertices, and adds the new angle to total_angle.

After it adds up the angles, the code returns True if total_angle is not zero. It also returns total_angle so the program can print it.

Calculating Angles

The get_radians function returns the angle between three points. It's unfortunate that the math library doesn't have a function to do this. Pygame's Vector2 class has an angle_to method that does this and I often use that class when I'm doing a lot of geometry, for example when writing games (see my book Build Your Own Python Action Arcade!), but for this I decided to stick with numbers.

The following code shows the get_radians function that the point_in_polygon method uses to find angles.

def get_radians(a, b, c): '''Return the angle ABC.''' # Return a value between PI and -PI. # Note that the value is the opposite of what you might # expect because Y coordinates increase downward. # Get the dot product. dot = dot_product(a, b, c) # Get the cross product. cross = cross_product_length(a, b, c) # Calculate the angle. return math.atan2(cross, dot)

This function uses the dot_product and cross_product_length functions to get the angle's dot product and length of the cross product.

If the two vectors are AB and BC, then the dot product is given by |AB| * |BC| * cos(theta) and the length of the cross product is given by |AB| * |BC| * sin(theta), where theta is the angle between the two vectors.

If you divide the length of the cross product by the length of the dot product, you get:

(|AB| * |BC| * sin(theta)) / (|AB| * |BC| * cos(theta)) = sin(theta)) / cos(theta)) = tan(theta)

The math.atan2 function takes opposite and adjacent side lengths for a right triangle and returns an angle with the corresponding tangent. That fact plus the previous equations means math.atan2(cross_product, dot_product) returns the angle between the two vectors.

Finding Dot Products

The previous section mentioned that the dot product is given by |AB| * |BC| * cos(theta), but there's also a direct way to calculate the dot product of two vectors. Simply multiply the vectors' corresponding X and Y coordinates and add them together. (For higher dimensional vectors, you multiply all of the corresponding coordinates.)

The following function calculates the dot product for the angle defined by three points.

def dot_product(a, b, c): '''Return the dot product ab · bc.''' # Note that ab · bc = |ab| * |bc| * cos(theta). # Get the vectors' coordinates. ba = (a[0] - b[0], a[1] - b[1]) bc = (c[0] - b[0], c[1] - b[1]) # Calculate the dot product. return ba[0] * bc[0] + ba[1] * bc[1]

This code subtracts the points to get the vectors pointing from point A to point B, and from point C to point B. It then multiples their corresponding X and Y coordinates and adds them together.

Finding Cross Products

As for the dot product, there's a simple way to calculate the cross product. It's more complicated so I'm not going to explain how to do it in higher dimensions, but it's easy enough for two-dimensional vectors.

The following cross_product_length function calculates the cross product.

def cross_product_length(a, b, c): '''Return the cross product AB × BC.''' # The cross product is a vector perpendicular to ab and bc having length # |ab| * |bc| * sin(theta) and direction given by the right-hand rule. # For two vectors in the X-Y plane, the result is a vector with X and Y # components 0 so the Z component gives the vector's length and direction. # Get the vectors' coordinates. ba = (a[0] - b[0], a[1] - b[1]) bc = (c[0] - b[0], c[1] - b[1]) # Calculate the Z coordinate of the cross product. return ba[0] * bc[1] - ba[1] * bc[0]

For more information on the cross product and how to calculate it, see Wikipedia or Wolfram MathWorld.

Now the get_radians function can use dot_product and cross_product_length to calculate angles.

Testing Overlapping Polygons

[This program draws points green if they're outside the polygon and red if they're inside it]

The discussion so far has assumed that the polygon is relatively simple like a rectangle or triangle. It still works if the polygon is complex (like the picture at the top of the post), but things are slightly different if the polygon is self-intersecting like the star on the right.

If the point is outside of the polygon, the same logic applies: clockwise and counterclockwise angles cancel each other out and you're left with a total angle of zero. If the point is inside the polygon, however, the sum of the angles may not be ±2π. Instead the total will be a multiple of ±2π that depends on the number of times the polygon winds around the point. In the picture on the right, the point under the crosshair is in the middle of the star so the polygon's points make 2 laps around the point and the sum is twice ±2π.

Polygons with more self-intersections can have even larger angle totals, but the total is always near zero if the point lies outside of the polygon. That means the point_in_polygon function shown here still works.

Conclusion

The point_in_polygon function lets you determine whether a point lies inside a polygon. It requires a few extra tools that find the angle dot product and the length of the cross product, which are somewhat confusing but easy to calculate. Once all the pieces are in place, the point_in_polygon function is easy to use.

Download the example to experiment with it and to see additional details.