Simple Point In Polygon
(convex) Tests

By Nathan Whitaker

 

Introduction

Some time ago I wanted a quick routine to detect whether a ray hit a polygon in 3 space. I was quite disappointed when I looked at the c.g.a faq and found that the recomended technique was to find the point of intersection of the ray in the polygons plane and then do a separate test to see if that point lay inside the polygon. With this in mind I present a simple one step technique harnassing the power of my favourite tool - the dot product.

 

2D Point In Poly Test

I have included the 2d test as it may be useful for such applications as selecting a polygon in a point and click environment etc.

This test involves generating a 2d normal from a plane which runs through each one of the polygons edges. Once we have this normal we take the dot product of the test point ( relative to any point on the edge - a vertex will do ) with the planes normal. If we build our normals so they point outwards from the center of the polygon then a positive value will mean that the test point lies outside of that edge and cant possibly be inside the polygon ( test over ). However if the value for all edges is negative the test point lies inside the polygon - simple! The only snag is that you will need two cases if you want to test polygons where the order of the polygons vertices can be in clockwise or anti-clockwise order. The difference between the two cases is simply the sign of the return value.

 

int point_in_poly( POINT *poly_verts, int num_verts, POINT *test_point )
{
	int nx, ny, loop;

//	Clockwise order.

	for ( loop = 0; loop < num_verts; loop++ )
	{
	//	generate a 2d normal ( no need to normalise ).
		nx = poly_verts[ ( loop + 1 ) % num_verts ].y - poly_verts[ loop ].y;
		ny = poly_verts[ loop ].x - poly_verts[ ( loop + 1 ) % num_verts ].x;

		x = test_point->x - poly_verts[ loop ].x;
		y = test_point->y - poly_verts[ loop ].y;

	//	Dot with edge normal to find side.
		if ( ( x * nx ) + ( y * ny ) > 0 )
			return 0;
	}

	return 1;
}

3D Ray Intersects Poly Test

The 3 space ray hits polygon test is based around the test above. If we generate a plane from the ray origin and a polygon edge we can perform a similiar test to the one above. In this example I assume that the ray extends to infinity from the ray origin along the direction formed from its start to end. Once again we need two cases dependant on the polygons orientation.

VECTOR CalcPlaneNormal( VECTOR *v0, VECTOR *v1, VECTOR *v2 )
{
        VECTOR normal;
        int x1, y1, z1, x2, y2, z2, x3, y3, z3;

        x1 = v0->vx;
        y1 = v0->vy;
        z1 = v0->vz;
        x2 = v1->vx;
        y2 = v1->vy;
        z2 = v1->vz;
        x3 = v2->vx;
        y3 = v2->vy;
        z3 = v2->vz;

        normal.vx = ( ( y1 - y2 ) * ( z1 - z3 ) ) - ( ( z1 - z2 ) * ( y1 - y3 ) );
        normal.vy = ( ( z1 - z2 ) * ( x1 - x3 ) ) - ( ( x1 - x2 ) * ( z1 - z3 ) );
        normal.vz = ( ( x1 - x2 ) * ( y1 - y3 ) ) - ( ( y1 - y2 ) * ( x1 - x3 ) );

        return normal;
}

// -----------------------------------------------------------------------------------------------------



int CheckRayPolyIntersection( int numverts, VECTOR *verts, VECTOR *raystart, VECTOR *rayend )
{
	VECTOR normal, v0, v1, relend;
	int loop, orientation;

//	Use z component of the plane normal to get orientation.

	orientation = ( ( verts[ 0 ].vx - verts[ 1 ].vx ) * ( verts[ 0 ].vy - verts[ 2 ].vy ) ) - 
                      ( ( verts[ 0 .vy 			- verts[ 1 ].vy ) * ( verts[ 0 ].vx -
                          verts[ 2 ].vx ) );

	if ( orientation > 0 )
	{
		for ( loop = 0; loop < numverts; loop++ )
		{
			v0.vx = verts[ loop ].vx;
			v0.vy = verts[ loop ].vy;
			v0.vz = verts[ loop ].vz;

			v1.vx = verts[ ( loop + 1 ) % num_verts ].vx;
			v1.vy = verts[ ( loop + 1 ) % num_verts ].vy;
			v1.vz = verts[ ( loop + 1 ) % num_verts ].vz;

			normal = CalcPlaneNormal( rstart, &v0, &v1 );

			relend.vx = rayend->vx - raystart->vx;
			relend.vy = rayend->vy - raystart->vy;
			relend.vz = rayend->vz - raystart->vz;

			if ( DotProduct( &relend, &normal ) < 0 )
				return 0;
		}

		return 1;
	}
	else
	{
		for ( loop = 0; loop < num_verts; loop++ )
		{
			v0.vx = verts[ loop ].vx;
			v0.vy = verts[ loop ].vy;
			v0.vz = verts[ loop ].vz;

			v1.vx = verts[ ( loop + 1 ) % num_verts ].vx;
			v1.vy = verts[ ( loop + 1 ) % num_verts ].vy;
			v1.vz = verts[ ( loop + 1 ) % num_verts ].vz;

			normal = CalcPlaneNormal( rstart, &v0, &v1 );

			relend.vx = rayend->vx - raystart->vx;
			relend.vy = rayend->vy - raystart->vy;
			relend.vz = rayend->vz - raystart->vz;

			if ( DotProduct( &relend, &normal ) > 0 )
				return 0;
		}

		return 1;
	}
}