Python: find if point lay on the border of a polygon

Question

I have a point-i and i wish to create a function to know if this point lies on the border of a polygon.

using:

def point_inside_polygon(x, y, poly):
    """Deciding if a point is inside (True, False otherwise) a polygon,
    where poly is a list of pairs (x,y) containing the coordinates
    of the polygon's vertices. The algorithm is called the 'Ray Casting Method'"""
    n = len(poly)
    inside = False
    p1x, p1y = poly[0]
    for i in range(n):
        p2x, p2y = poly[i % n]
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xinters = (y-p1y) * (p2x-p1x) / (p2y-p1y) + p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x, p1y = p2x, p2y
    return inside

I am able to know only if the points lies within the polygon.

poly = [(0,0), (2,0), (2,2), (0,2)]
point_inside_polygon(1,1, poly)
True
point_inside_polygon(0,0, poly)
false
point_inside_polygon(2,0, poly)
False
point_inside_polygon(2,2, poly)
True
point_inside_polygon(0,2, poly)
True

How can I write a function to find if a point lay on the border of a polygon instead?

Answer 1

Everyone is overcomplicating things. Here is a short point on polygon, assuming you have a distance function and a small EPSILON.

def pointOnPolygon(point, poly):
    for i in range(len(poly)):
        a, b = poly[i - 1], poly[i]
        if abs(dist(a, point) + dist(b, point) - dist(a, b)) < EPSILON:
            return true
    return false

Answer 2

Found the solution to this at: http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html

Here the code:

# Improved point in polygon test which includes edge
# and vertex points

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly[i]
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "IN"

   n = len(poly)
   inside = False

   p1x,p1y = poly[0]
   for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
         if y <= max(p1y,p2y):
            if x <= max(p1x,p2x):
               if p1y != p2y:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "IN"
   else: return "OUT"

# Test a vertex for inclusion
polygon = [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat= -33.416032
lon= -70.593016

print point_in_poly(lat, lon, polygon)

# test a boundary point for inclusion
poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)]
x = 3
y = 1
print point_in_poly(x, y, poly2)

Answer 3

For each pair of adjacent vertices A,B:

1. construct a vector from A to B, call it p

2. now construct a vector from A to your test point X call it q

3. the dot product formula for a pair of vectors is p.q = |p||q|cosC where C is the angle between the vectors.

4. so if p.q/|p||q| == 1 then the points AX and AB are co-linear. Working on a computer, you will want 1 - p.q/|p||q| < some_small_value depending on how accurate you want to be.

5. also need to check that |q| < |p| (ie X is closer to A than B)

if 4&5 are true your point is on the border.

Edit

The other way I think I've seen this done is to take your test point X, and construct a line through X perpendicular to the line between A and B. Find where this line and the line A->B cross. Work out the distance from X to this crossing point, if that is sufficiently small you consider the point as being on the line.

Edit -- fun little exercise!

Posted some code that was wrong earlier due to me misremembering some maths. Had a play in Pythonista on the train home and came up with this which seems to basically work. Have left the maths proof out as editing posts on iPad is painful!

Not much testing done, no testing for division by zero etc, caveat user.

    # we determine the point of intersection X between
    # the line between A and B and a line through T
    # that is perpendicular to the line AB (can't draw perpendicular
    # in ascii, you'll have to imagine that angle between AB and XT is 90
    # degrees.
    #
    #       B
    #      /
    #.    X  
    #    / \
    #   /   T
    #  A
    # once we know X we can work out the closest the line AB
    # comes to T, if that distance is 0 (or small enough)
    # we can consider T to be on the line
    import math


    # work out where the line through test point t
    # that is perpendicular to ab crosses ab
    #
    # inputs must be 2-tuples or 2-element lists of floats (x,y)
    # returns (x,y) of point of intersection
    def intersection_of_perpendicular(a,b,t):

    if a[0] == b[0]:
            return (a[0],t[1])

    if a[1] == b[1]:
            return (t[0],a[1])

    m = (a[1] - b[1])/(a[0] - b[0]) #slope of ab

    x_inter = (t[1] - a[1] + m*a[0] + (1/m)*t[0])*m/(m**2 + 1)
    y_inter = m*(x_inter - a[0]) + a[1]
    y_inter2 = -(1/m)*(x_inter - t[0]) + t[1]

    #print '...computed ',m,(x_inter, y_inter), y_inter2
    return (x_inter, y_inter)

    # basic Pythagorean formula for distance between two points
    def distance(a,b):
        return math.sqrt( (a[0]-b[0])**2 + (a[1]-b[1])**2 )

    # check if a point is within the box defined by a,b at
    # diagonally opposite corners
    def point_in_box(a,b,t):
        xmin = min(a[0],b[0])
        xmax = max(a[0],b[0])
        ymin = min(a[1],b[1])
        ymax = max(a[1],b[1])

        x_in_bounds = True
        if xmax != xmin:
            x_in_bounds = xmin <= t[0] <= xmax
        y_in_bounds = True
        if ymax != ymin:
            y_in_bounds = ymin <= t[1] <= ymax
        return x_in_bounds and y_in_bounds

    # determine if point t is within 'tolerance' distance
    # of the line between a and b
    # returns Boolean
    def is_on_line_between(a,b,t,tolerance=0.01):
        intersect = intersection_of_perpendicular(a,b,t)
        dist = distance(intersect, t)
        in_bounds = point_in_box(a,b,t)
        return in_bounds and (dist < tolerance)


a = (0,0)
b = (2,2)
t = (0,2)

p = intersection_of_perpendicular(a,b,t)
bounded = point_in_box(a,b,t)
print 'd ',distance(p,t), ' p ',p, bounded

a = (0,2)
b = (2,2)
t = (1,3)

p = intersection_of_perpendicular(a,b,t)
bounded = point_in_box(a,b,t)
print 'd ',distance(p,t),' p ',p, bounded

a = (0.0,2.0)
b = (2.0,7.0)
t = (1.7,6.5)

p = intersection_of_perpendicular(a,b,t)
bounded = point_in_box(a,b,t)
on = is_on_line_between(a,b,t,0.2)
print 'd ',distance(p,t),' p ',p, bounded,on 

Answer 4

For a convex polygon, you can solve the problem in O(log n) time by ordering vertices clock-wise and storing for each vertex the angle between the vertex and a point c in the interior of the polygon. Then for a query point x, you get the angle from c to x and binary search to find the unique pair of adjacent vertices (v1,v2) such that the angle to x is between the angles to v1 and to v2. Then x is either on the edge (v1,v2), or x is not on the boundary.

If you have a more complicated polygon, then you can try decomposing the polygon into a union of convex polygons by adding some interior edges (e.g. first triangulate and then remove edges to get bigger convex polygons); if the number of convex polygons you get is small (say k), then you can test each convex polygon to see if a point is on an edge, and the overall running time is O(k lg n) where n is the total number of vertices in your polygon.

Alternatively, if you are not worried about using extra space, and you really want to quickly determine if you're on an edge, then you can break up each edge into equally spaced segments by adding extra "vertices" along each edge; it's hard to say how many is enough (sounds like an interesting math problem), but clearly if you add enough extra vertices along each edge then you can tell which edge a point must lie on simply by finding the nearest neighbor out of your set of vertices (original vertices and the ones you added), and then just test the one or two edges that the nearest neighbor vertex lies on. You can very quickly find k-nearest neighbors in 2-d if you use a 2-dimensional kd-tree (you build the tree as a pre-processing step, and then the tree supports fast k-nearest neighbor queries), and the kd-tree tree only uses linear space.

Answer 5

It might help to break down the problem into three steps:

1. Write a function that can determine if a point is on a line segment.
2. Compute all of the line segments that make up the border of the polygon.
3. The point is on the border if it is on any of the line segments.

Here's some python code, assuming you've written or found a suitable candidate for isPointOnLineSegmentBetweenPoints:

def pointOnPolygon(point, polygonVertices):
    n = len(polygonVertices)
    for i in range(n):
        p1 = polygonVertices[i]
        p2 = polygonVertices[-n+i+1]
        if isPointOnLineSegmentBetweenPoints(point, p1, p2):
            return true
    return false

Answer 6

I haven't tested this, but the general idea is:

def pointOnBorder(x, y, poly):
    n = len(poly)
    for(i in range(n)):
        p1x, p1y = poly[i]
        p2x, p2y = poly[(i + 1) % n]
        v1x = p2x - p1x
        v1y = p2y - p1y #vector for the edge between p1 and p2
        v2x = x - p1x
        v2y = y - p1y #vector from p1 to the point in question
        if(v1x * v2y - v1y * v2x == 0): #if vectors are parallel 
            if(v2x / v1x > 0): #if vectors are pointing in the same direction
                if(v1x * v1x + v1y * v1y >= v2x * v2x + v2y * v2y): #if v2 is shorter than v1
                    return true
    return false