Check if a point is within a polygon – general case

Previously we had shown how to check if a two dimensional point is within a convex polygon. Now let’s look at the general case, where the polygon can be either convex or non-convex. To do so, we’ll use the ray casting algorithm. Briefly, this algorithm creates a ray or line segment from a point outside of the polygon to the point in question. Counting up the number of intersections, if the count is even (or zero), then the point lies outside of the polygon. Alternatively, if the count is odd, then the point is within the polygon.

So let’s create the general function pointWithinPolygon.m, which you can download and try yourself, to check if a given point is within a polygon defined by x and y coordinates px and py. This function will utilize the previously described checkSegmentIntersection and checkPointOnSegment functions. Once again we will define the polygon in either the clockwise or counterclockwise directions without repeating the first vertex. For the polygon, let’s define the following non-convex polygon as follows:

px = [0, 1,-1,-1,2,2,3,4];
py = [0, -1,-1,1,2,4,4,1];

First, let’s initialize our function and do some error checking to ensure proper definitions for the polygon and point of interest:

function insidePoly = pointWithinPolygon(px, py, point, plotResults)

if numel(px) ~= numel(py) || numel(px) < 3
    error('Incorrectly defined Polygon');
end
if numel(point) ~= 2
    error('Incorrectly defined point');
end

Now, to generate a count of how many segment intersections occur for a given ray from outside the polygon to the point in question, we first need to generate an outside point and the ray. To ensure that the outsidePoint is not within the polygon, we will utilize the minimum values minus one along all x and y of the polygon. The ray will be defined as an array from the defined outsidePoint to the point in question.

%% Generate ray from outside polygon to the point in question.
outsidePoint = [min(px) - 1, min(py) - 1];
ray = [outsidePoint; point];

The count for segment intersection will performed by looping through each edge of the polygon and checking if the given segment and ray will intersect. The intersection check will once again use the previously described checkSegmentIntersection code. A true value indicates that the ray intersects a segment and will lead to an increment in the count.

%% Check the intersection of each segment with the ray 
count = 0;
for k = 1:numel(px)
    if k < numel(px)
        doesIntersect = checkSegmentIntersection(ray, [px(k) py(k); px(k+1) py(k+1)], false);
    else
        doesIntersect = checkSegmentIntersection(ray, [px(k) py(k); px(1) py(1)], false);
    end
    if doesIntersect
        count = count + 1;
    end
end

This code is almost perfect… Unfortunately something goes horribly wrong when the ray passes through vertices of the polygon. For example, if we want to check the point (1,1) and select the outside point from (-2,-2), we will pass through the vertices (-1,-1) and (0,0), resulting in a count of 4, and incorrectly predicting that the point is outside the polygon:

To fix this problem, we’ll check if any vertices lie on the ray (using the code provided previously in checkPointOnSegment), and if they do, rotate the outside point until no vertices intersect the ray. This part of the code will have to be performed prior to the previously posted segment intersection check.

% ensure that ray does not intersect with polygon vertices and modify the
% outside point if it does so.
checkVert = true;
while checkVert
    for k = 1:numel(px)
        onRay = checkPointOnSegment(ray, [px(k), py(k)], false);
        if onRay
            outsidePoint = outsidePoint - rand(1,2);
            ray = [outsidePoint; point];
            break
        end
        if k == numel(px)
            checkVert = false; 
        end
    end
end

Now, our code should be robust with example results provided below. Please give the pointWithinPolygon.m code a try and let us know if you have any comments!

Check if a point is on a line segment

As we continue to look at different computational geometry problems, one simple, yet important test is to see whether a given point lies on a line segment. Similar to our previous posts, we will be utilizing the cross product and dot product to check for the proper conditions. Hopefully, these computations have become second nature now. Again, I have uploaded the script checkPointOnSegment.m if you want to test on your own.

We will first define our segments and point of interest. Our segment will consist of a 2×2 array, with each row specifying each endpoint A and B, and a point C will be defined as a 2 element array. For example, given the following, we have a line from (-1,4) to (3,-2) and want to check if the point (1,1) lies on the segment.

>> segment = [-1 4; 3 -2];
>> C = [1 1];

Let’s now define the line segment as a vector from point A (-1,4) to point B(3,-2) and the vector AC:

A = segment(1,:);
B = segment(2,:);
% form vectors for the line segment (AB) and the point to one endpoint of
% segment
AB = B - A;
AC = C - A;

As we have shown previously, collinearity can be checked through the cross product. If the cross product of AB and AC is 0, then the segments have a 0 or 180 degree angle between them and are therefore collinear. We can do this simple check for whether a point is on a line, given our definition of cross. Note: if you use Matlab’s built-in cross, you’ll have to append 0’s to all of your vectors, but you will still obtain the same results.

% if not collinear then return false as point cannot be on segment
if cross(AB, AC) == 0
    checkPt = true;
end
%% Cross product returning z value
function z = cross(a, b)
    z = a(1)*b(2) - a(2)*b(1);

Alas, a line segment is different than a line, in that it has endpoints, while a line continues to infinity in both directions. So we need to perform some additional checks for whether the point is beyond the limits of the line segment. While we could do some simple > or < checks, let's use the dot product to see whether the point lies on the line segment or is found at the endpoints. So our complete code will be:

% if not collinear then return false as point cannot be on segment
if cross(AB, AC) == 0
    % calculate the dotproduct of (AB, AC) and (AB, AB) to see point is now
    % on the segment
    dotAB = dot(AB, AB);
    dotAC = dot(AB, AC);
    % on end points of segment
    if dotAC == 0 || dotAC == dotAB
        onEnd = true;
        checkPt = true;
    % on segment
    elseif dotAC > 0 && dotAC < dotAB
        checkPt = true;
    end
end

Let’s look at some examples. Please feel free to try on your own and let us know how it works for you.

[checkPt, onEnd] = checkPointOnSegment(segment, C, true)
checkPt =
  logical
   1

onEnd =
  logical
   0


[checkPt, onEnd] = checkPointOnSegment(segment, [7,-5], true)
checkPt =
  logical
   0

onEnd =
  logical
   0