which side of a polyline a point is

Question

I'm trying to write a code that would allow to identify which side of a random polyline (not closed one) a random point is. Polyline is considered infinite with the first and last segments extended to infinity. If they are crossed polyline should be considered as a polygon (I haven't elaborate the code for this case yet). The logic is following:

1. Define the vertex of a polyline, the distance to which from the point under consideration is minimum. Variable minimum2p is the distance to that vertex, I2p is the index of the vertex.

2. Define the segment of a polyline, the distance to which from the point under consideration is minimum. Only those segments count (variable count_s) which are intersected by a perpendicular from the point under consideration. Variable minimum2s is the minimum distance to the segment; I2p is the index of the first vertex of that segment; flag is the boolean variable, which stores the information on which segment is intersected by the mentioned perpendicular.

3. Next it's just a matter to choose proper segment to compare with using, for example, ideas from the links link-1, link-2 or link-3. I tried the approach from here, but it doesn't work for many special cases. I use the best answer there for internal points of a polyline though. So, my approach is the following:

4. First, check if it's the first or the last vertex of a polyline. If this is the case, then the chosen segment is either first or last one correspondingly, but only if there is no other segment closer than the first or the last ones. If there is another segment then I chose that segment.

5. Next, if step 4 is not the case, then I check internal vertex of a polyline. If there is also a segment close nearby then I compare indexes I2p and I2s, if the last exists. If they coincide then there are no ambiguities in choosing the proper segment to compare with. If they are different then preference goes to the closest segment rather than to the closest vertex.

6. Finally, if there is no segment nearby (in the sense of perpendicular from the point crossing the segment), then for the internal vertex I apply the idea from the best answer here.

Here are some results for different polylines, which are defined by the X and Y coordinates of their vertexes, stored in 'polylineX' and 'polylineY' correspondingly (red color is for the 'left' position, grey color is for the 'right' position, black color is the position on the polyline, blue line represents a polyline).

As you can notice, for relatively smooth polylines the code works. However, for sharper ones, or in some ways complicated ones, the code doesn't work properly. What do I miss in my code? What condition should be added to take some cases into account?

The code is the following:

clear all
close all 
clc
clf
polylineX = [0 1 2 3 4 5 6 7 8];
polylineY = [-10 20 -13 18 -17 16 -21 23 -25];
hold on
title(['polylineX=[',num2str(polylineX),'], polylineY=[',num2str(polylineY),']'])
chosen = 0;
span = 60;

for ii = 10:70
for jj = 30:60
    ii
    jj
    position = -2;
    point = [(jj-round(span/2))/1 (ii-round(span/2))/1];

    axis equal
    plot(polylineX,polylineY,'.-','MarkerSize',1,'LineWidth',1);

    distance2p = zeros(1,length(polylineX)); % distances from the point to the points (2p) of the polyline
    distance2s = zeros(1,length(polylineX)-1); % distances from the point to the segments (2s) of the polyline
    flag = zeros(1,length(polylineX)-1);
    count_s = 0; % counter of segments, which are intersected by the normal pointing from the 'point'
    k = 0;

    for i = 1:length(polylineX)-1
        pos = sign((polylineX(i+1) - polylineX(i)) * (point(2) - polylineY(i)) -...
                   (polylineY(i+1) - polylineY(i)) * (point(1) - polylineX(i)));

        % computing the distances from the 'point' to all segments and mark if
        % the distance vectors intersect the segments
        [flag(i),distance2s(i)] = distanceToLine([polylineX(i) polylineX(i+1)],[polylineY(i) polylineY(i+1)],[point(1) point(2)]);

        if flag(i)
            if k == 0
                minimum2s = distance2s(i);
                I2s = i;
            end;
            k = 1;
            count_s = count_s + 1; % count segments, which are intersected by the normal pointing from the 'point'
            if distance2s(i) < minimum2s
                I2s = i;
                minimum2s = distance2s(i);
            end;
        end;
    end;


    % first compute the distances between the 'point' under consideration and the
    % points of the given polyline
    for i  = 1:length(polylineX)
        distance2p(i) = sqrt((point(1)-polylineX(i))^2+(point(2)-polylineY(i))^2);
    end;
    [minimum2p,I2p] = min(distance2p);
    clear k pos i

    % now we need to choose which segment of the polyline to compare our 'point' with. These
    % segments are either adjacent to that point of the polyline, which is the closest
    % to the 'point' of interest, or the closest to the 'point' segment, which
    % has an intersection with the normale pointing from the 'point'.

    if I2p == 1 % if the 'point' is near the start of polyline
        if exist('minimum2s','var')
            if I2p == I2s
                chosen = I2p;
            else
                chosen = I2s;
            end;
        else
            chosen = I2p;
        end;

    elseif I2p == length(polylineX) % if the 'point' is near the end of polyline
        if exist('minimum2s','var')
            if I2s == I2p-1
                chosen = I2p - 1;
            else
                chosen = I2s;
            end;
        else
            chosen = I2p - 1;
        end;
    else
        if exist('minimum2s','var')
            if I2p == I2s
                chosen = I2p;

            else
                chosen = I2s;
            end;
        else
                pos1 =  sign((polylineX(I2p) - polylineX(I2p-1)) * (point(2) - polylineY(I2p-1)) -...
                 (polylineY(I2p) - polylineY(I2p-1)) * (point(1) - polylineX(I2p-1)));
                % position of the second segment relative to the first segment
                pos2 =  sign((polylineX(I2p) - polylineX(I2p-1)) * (polylineY(I2p+1) - polylineY(I2p-1)) -...
                             (polylineY(I2p) - polylineY(I2p-1)) * (polylineX(I2p+1) - polylineX(I2p-1)));
                if (pos1 == 1 && pos2 == 1) || (pos1 == -1 && pos2 == -1)
                    chosen = I2p;
                elseif pos1 == 0 || pos2 == 0
                    chosen = I2p;
                else
                    chosen = I2p - 1;
                end;
        end;
    end;

    position = sign((polylineX(chosen+1) - polylineX(chosen)) * (point(2) - polylineY(chosen)) -...
                    (polylineY(chosen+1) - polylineY(chosen)) * (point(1) - polylineX(chosen)));

    if position == 1
        plot(point(1),point(2),'r.','MarkerSize',5)
    elseif position == -1;
        plot(point(1),point(2),'.','Color',[0.9 0.9 0.9],'MarkerSize',5) % gray color
    elseif position == 0
        plot(point(1),point(2),'k.','MarkerSize',5)
    elseif position == -2
        plot(point(1),point(2),'g.','MarkerSize',5)
    end;

    pause(0.00000001)
    clear chosen  count_s distance2p distance 2s flag I2p I2s minimum2p minimum2s point pos1 pos2 position

end;
end;

Answer 1

Naive idea I have is to integrate angle from point to line. Integrations starts from one side at infinity then over all points to another side of infinity. That is just bunch of atan2 functions.

Side of curve is determinated from sign of integration. This should work even if curve is overlapping.

In python:

from math import atan2,pi
#import matplotlib.pyplot as plt

# difference of two angles should be always -pi < x < pi
def fixed_angle(a):
    if a > pi:
       return a - 2*pi
    elif a < (-1*pi):
       return a + 2*pi
    assert(-1*pi < a < pi) # just check, to be sure
    return a

# polyline
xs = [0, 1, 2, 3, 4, 5, 6, 7, 8];
ys = [-10, 20, -13, 18, -17, 16, -21, 23, -25];

# start point
x = 4
y = 0

#from first two points
angle_start = atan2(ys[0]-ys[1],xs[0]-xs[1])
#last angle is angle of last section
angle_end = atan2(ys[-1]-ys[-2],xs[-1]-xs[-2]) 

integ = 0
prev =  angle_start
for i in range(len(xs)):
    a = atan2(ys[i]-y,xs[i]-x)
    integ += fixed_angle(a-prev)
    prev = a

integ += fixed_angle(angle_end - prev)

if integ > 0:
    print("point is left")
else:
    print("point is right")

#plt.plot(xs,ys)
#plt.show()

Answer 2

How about an alternative solution, based on something like the notion of "winding number"? This would involve computing the total angle subtended by all segments of your boundary polyline as viewed from your test point.

If you imagine that this polyline is extended by a semicircle at infinity to make a closed contour, then this sum of angles will be zero for a point outside the contour but 2*Pi for a point inside. Obviously, whether the semicircle is in the left or right half-plane will determine which side of the line is "inside" or "outside". If you exclude the semicircle itself from the angle-sum, because it contributes +/-Pi, then a point on one side of the line should have an angle-sum of +Pi, while the other side should give -Pi. Whether it's +Pi or -Pi depends on the direction defined by the sequence of vertices in your polyline.

One subtlety with this approach is computing the (signed) angle subtended by each segment of the polyline. I believe this could be done by looking at the two vectors that join the test point to the two ends of the given segment. By taking the dot-product of those vectors (and after dividing by their lengths), you can compute the cosine of the angle subtended by that segment. Taking the inverse-cosine will be ambiguous between positive of negative angles. However, by taking the cross-product of the two vectors (treating their z-components as zero) you can compute the sine of the angle, whose sign will allow you to compute the correct branch of the inverse-cosine.

Answer 3

An approach that doesn't calculate angles is to consider the signed areas of triangles made by the test point and successive polyline segments.

The signed area of a triangle can be understood (by a convention you adopt) as one vertex being to the 'left' or 'right' of the opposite line. It's calculated by the formula 0.5 * (-x2*y1 + x3*y1 + x1*y2 - x3*y2 -x1*y3 + x2*y3) for three vertices (xi, yi). Crucially, the order you traverse the vertices affects the sign.

Suppse your polyline is a list [v1, v2, ...., vn]. For each ordered pair of adjacent vertices (v_i, v_(i+1)), compute the signed area of triangle (v_i, your-test-point, v_(i+1)). It seems that the one to care about for 'side' decisions is the triangle with the smallest absolute area: that is, where the point is 'closest' to the polyline. Decide by your convention whether that triangle is left- or right-side with respect to your polyline by considering its signed area.

*edit - a zero-area triangle means your test point is colinear with a polyline segment; you'd need to test whether it was actually on the segment separately.

Answer 4

my solution is self-explained in the figure.

1. project the point, e.g. P₁ to the polyline [O₁,O₂,O₃,O₄] and the projection point is P^'₁ which is at the line O₁O₂.
2. If O₁P^'₁P₁ and P₁P^'₁O₂ is both counter clockwise oriented, then the point P₁ is at the left side of the polyline.
3. If the point (e.g. P₂) at the right side of the polyline, then O₃P^'₂P₂ and P₂P^'₂O₄ are both clockwise oriented.