Point in Polygon Test

Given a polygon and a point, how do you tell if the point is inside or outside the polygon?

Idea

In 2D let’s think of the polygon as dividing the whole space into two regions. The finite region inside the polygon and the infinite region outside. Imagine yourself standing at the point which we wish to classify as being inside or outside. Next, pick a direction and start walking in a straight line. As you walk, keep track of the number of times you cross an edge of the polygon. After a long while when you are certain you have left the polygon behind, maybe at infinity, take a break and consider the number of times you crossed its edges.

Every edge crossing implies that you have either left or entered the polygon so it toggles your current state of being inside or outside. An odd number of crossings means that your current state is opposite of what it was when you started while an even number of crossings means it is the same. So as you take your break at infinity far outside of the polygon you know that if the number of crossings was odd you started on the inside and if it was even you started outside.

Implementation

Now let’s convert the idea into some code. Here is what a polygon looks like:

struct Polygon
{
    vec2 *Vertices;
    int   NumVertices;
};

A polygon is a list of vertices. The edges are implied by the order of the vertices with each pair of consecutive vertices forming an edge. In this case the algorithm translates to code fairly naturally:

// This function returns true if the ray with origin o and direction d intersects a line segment
// with endpoints a and b.
bool RayLineSegmentIntersection( const vec2 &o, const vec2 &d, const vec2 &a, const vec2 &b );

bool IsInside( const vec2 &point, const Polygon &polygon )
{
    int numCrossings = 0;

    // Loop over all edges in the polygon.
    for( int i = 0; i < polygon.NumVertices; ++i )
    {
        int j = ( i + 1 ) % polygon.NumVertices;

        // Check if a ray in the positive x direction crosses the current edge.
        if( RayLineSegmentIntersection( point, vec2( 1, 0 ), polygon.Vertices[ i ], polygon.Vertices[ j ] ) )
            ++numCrossings;
    }

    // If the number of crossings is odd, the point was inside the polygon.
    return ( numCrossings % 2 ) == 1;
}

The implementation above does exactly what we described. It shoots a ray from the point in question down the x axis(the direction doesn’t matter since a polygon will completely surround us if we’re inside) and counts how many edges it intersects. Then depending on the count the point is classified.

Ray-Line Segment Intersection Test

Here is a way to test a ray against a line segment. To get this result just set the ray \mathbf{x_1}(t_1)=\mathbf{o}+\mathbf{d}t_1 equal to the line segment \mathbf{x_2}(t_2)=\mathbf{a} + (\mathbf{b} - \mathbf{a})t_2 and solve for t_1, t_2. The ray will hit the line segment if t_1 \geq 0(line segment is in front of the ray) and 0 \leq t_2 \leq 1(the ray crosses the line segment between \mathbf{a} and \mathbf{b}).

bool RayLineSegmentIntersection( const vec2 &o, const vec2 &d, const vec2 &a, const vec2 &b )
{
    vec2 ortho( -d.y, d.x );
    vec2 aToO( o - a );
    vec2 aToB( b - a );

    float denom = dot( aToB, ortho );

    // Here would be a good time to see if denom is zero in which case the line segment and
    // the ray are parallel.

    // The length of this cross product can also be written as abs( aToB.x * aToO.y - aToO.x * aToB.y ).
    float t1 = length( cross( aToB, aToO ) ) / denom;
    float t2 = dot( aToO, ortho ) / denom;

    return t2 >= 0 && t2 <= 1 && t1 >= 0;
}

That is all there is to it. In a future post we’re going to revisit and solve this problem using some tricks!

Advertisements
Point in Polygon Test

14 thoughts on “Point in Polygon Test

    1. Scala says:

      And I’m curious why you hard code the second point on the Ray: vec2( 1, 0 )
      if( RayLineSegmentIntersection( point, vec2( 1, 0 ), polygon.Vertices[ i ], polygon.Vertices[ j ] ) )

      1. That’s the direction of the ray. Since we only care if the ray leaves the polygon the direction doesn’t matter. I just happened to pick the unit x vector.

  1. eliw00d says:

    I am a little confused about “length( cross( aToB, aToO ) )”. Is there a length function you are using, or is this pseudo code? Your comment just gives the definition for cross product, so it is a little unclear as to what should be done in code there.

    Also, would the intersection point be “new Vector2(o.x + (d.x * t1), o.y + (d.y * t1))”?

    1. The code uses glm as the math library. The length function just gives you the magnitude of the vector. You’re right to be confused! I guess I was thinking in 3D when writing this bit and used the library functions because it looked neater. You can write float t1 = abs(aToB.x * aToO.y - aToO.x * aToB.y) / denom. As for the intersection point, you’re right, that will give you the position along the ray.

      1. eliw00d says:

        Awesome, thanks! One other question: I seem to be getting intersections that would only happen if intersecting rays.

        Any ideas why that might be happening with this algorithm?

      2. Hey! This whole page has been super helpful to me, but this line kinda threw me for a loop. For reasons I don’t fully understand, this implementation didn’t work for me until I removed the “abs”. Any idea why that would be?

    1. That makes sense because that step of the algorithm uses the length of the cross product and the length is always positive. I’ll update the comment as it is misleading! Thanks!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s