Create non-intersecting polygon passing through all given points

Question

Suppose I have an array of points in random order, and I need to find a polygon (by sorting them, such that every adjacent pair represents a side) which passes through all of the points, and its sides are non-intersecting of course.

I tried to do it by selecting a point, and adding all points to the final array which are below it, sorted left to right. Then, adding all points which are above it, sorted right to left.

I've been told that I can add an additional point and sort naturally to avoid self-intersections.. I am unable to figure out that though. What's a simple way to do this?

Answer 1

Our strategy is to make a plan where we are sure that the polygon includes all points, and that we can find an order to connect them where none of the lines intersect.

Algorithm:

1. Find the leftmost points p
2. Find the rightmost point q
3. Partition the points into A, the set of points below pq, and B, the set of points above pq [you can use the left turn test on (p,q,?) to determine if a point is above the line].
4. Sort A by x-coordinate (increasing)
5. Sort B by x-coordinate (decreasing).
6. Return the polygon defined by p, the points in A, in order, q, the points of B in order.

Runtime:

Steps 1,2,3 take O(n) time.
Steps 4,5 take O(nlogn) time.
Step 6 take O(n) time.
Total runtime is O(nlogn).

Correctness:

By construction, all points besides p,q are in set A or set B. Our output polygon from line 6 therefore outputs a polygon with all the points. We now need to argue that none of the line segments in our output polygon intersect each other.

Consider each segment in the output polygon. The first edge from p to the first point in A can't intersect any segment (because there is no segment yet). As we proceed in order by x-coordinate through the points in A, from each point, the next segment is going to the right, and all previous segments are to the left. Thus, as we go from p, through all the points of A, to point q, we will have no intersections.

The same is true as we go from q back through the points of B. These segments cannot intersect each other because they proceed from right to left. These segments also cannot intersect anything in A because all points in A are below line pq, and all points in B are above this line.

Thus, no segments intersect each other and we have a simple polygon.

Source: Broken link

Answer 2

Warning! Sometimes polygons intersect, I don't know why. This could be my implementation problem. See comments for intersection examples.

Check this answer before using my code: https://stackoverflow.com/a/64459159/2829863

Here is python 3.6 code (libraries required: matplotlib, numpy) based on bdean20's answer.

Pictures description:

• Top left - predefined polygon, other polygons are randomly generated.
• Dotted line - connects green (leftmost) and red (rightmost) polygon's points.
• Black dots are lays on the dotted line.
• Orange dots lays below dotted line.
• Blue dots lays above dotted line.

=========

import random
from operator import itemgetter
import numpy
import matplotlib
import matplotlib.pyplot

class Create_random_polygon:
    
    def __init__(self, array, min_rand_coord = None, max_rand_coord = None, points_num = None):        
        self.array = array
        self.min_rand_coord = min_rand_coord 
        self.max_rand_coord = max_rand_coord
        self.points_num = points_num

    def generate_random_points(self):
        random_coords_list = []
        for x in range(self.points_num):
            coords_tuple = (random.randint(self.min_rand_coord, self.max_rand_coord),
                            random.randint(self.min_rand_coord, self.max_rand_coord))
            random_coords_list.append(coords_tuple)
        self.array = random_coords_list
        return random_coords_list
    
    def close_line_to_polygon(self):
        a = self.array[0]
        b = self.array[len(self.array)-1]
        if a == b:
            pass
        else:
            self.array.append(a)    

    def find_leftmost_point(self):
        leftmost_point = None
        leftmost_x = None
        for point in self.array:
            x = point[0]
            if leftmost_x == None or x < leftmost_x:
                leftmost_x = x
                leftmost_point = point
        return leftmost_point

    def find_rightmost_point(self):
        rightmost_point = None
        rightmost_x = None
        for point in self.array:
            x = point[0]
            if rightmost_x == None or x > rightmost_x:
                rightmost_x = x
                rightmost_point = point
        return rightmost_point

    def is_point_above_the_line(self, point, line_points):
        """return 1 if point is above the line
           return -1 if point is below the line
           return  0 if point is lays on the line"""
        px, py = point
        P1, P2 = line_points
        P1x, P1y = P1[0], P1[1]
        P2x, P2y = P2[0], P2[1]
        array = numpy.array([
            [P1x - px, P1y - py],
            [P2x - px, P2y - py],
            ])
        det = numpy.linalg.det(array)
        sign = numpy.sign(det)
        return sign
    
    def sort_array_into_A_B_C(self, line_points):
        [(x_lm, y_lm), (x_rm, y_rm)] = line_points
        A_array, B_array, C_array = [], [], []
        for point in self.array:
            x, y = point
            sing = self.is_point_above_the_line( (x, y), line_points)
            if sing == 0:
                C_array.append(point)
            elif sing == -1:
                A_array.append(point)
            elif sing == 1:
                B_array.append(point)
        return A_array, B_array, C_array

    def sort_and_merge_A_B_C_arrays(self, A_array, B_array, C_array):
        A_C_array = [*A_array, *C_array]
        A_C_array.sort(key=itemgetter(0))
        B_array.sort(key=itemgetter(0), reverse=True)        
        merged_arrays = [*A_C_array, *B_array]
        self.array = merged_arrays

    def show_image(self, array, line_points, A_array, B_array, C_array):
        [(x_lm, y_lm), (x_rm, y_rm)] = line_points        
        x = [x[0] for x in array]
        y = [y[1] for y in array]
        Ax = [x[0] for x in A_array]
        Ay = [y[1] for y in A_array]
        Bx = [x[0] for x in B_array]
        By = [y[1] for y in B_array]
        Cx = [x[0] for x in C_array]
        Cy = [y[1] for y in C_array]          
        matplotlib.pyplot.plot(Ax, Ay, 'o', c='orange') # below the line
        matplotlib.pyplot.plot(Bx, By, 'o', c='blue') # above the line
        matplotlib.pyplot.plot(Cx, Cy, 'o', c='black') # on the line
        matplotlib.pyplot.plot(x_lm, y_lm, 'o', c='green') # leftmost point
        matplotlib.pyplot.plot(x_rm, y_rm, 'o', c='red') # rightmost point
        x_plot = matplotlib.pyplot.plot([x_lm, x_rm], [y_lm, y_rm], linestyle=':', color='black', linewidth=0.5) # polygon's division line
        x_plot = matplotlib.pyplot.plot(x, y, color='black', linewidth=1) # connect points by line in order of apperiance        
        matplotlib.pyplot.show()

    def main(self, plot = False):
        'First output is random polygon coordinates array (other stuff for ploting)'
        print(self.array)
        if self.array == None:
            if not all(
                [isinstance(min_rand_coord, int),
                 isinstance(max_rand_coord, int),
                 isinstance(points_num, int),]
                ):
                print('Error! Values must be "integer" type:', 'min_rand_coord =',min_rand_coord, ', max_rand_coord =',max_rand_coord, ', points_num =',points_num)
            else:                
                self.array = self.generate_random_points()            

        print(self.array)
        x_lm, y_lm = self.find_leftmost_point()
        x_rm, y_rm = self.find_rightmost_point()
        line_points = [(x_lm, y_lm), (x_rm, y_rm)]

        A_array, B_array, C_array = self.sort_array_into_A_B_C(line_points)
        self.sort_and_merge_A_B_C_arrays(A_array, B_array, C_array)
        self.close_line_to_polygon()
        if plot:
            self.show_image(self.array, line_points, A_array, B_array, C_array)
        return self.array

if __name__ == "__main__":
    # predefined polygon
    array = [ 
        (0, 0),
        (2, 2),
        (4, 4),
        (5, 5),
        (0, 5),        
        (1, 4),
        (4, 2),
        (3, 3),
        (2, 1),
        (5, 0),
        ]    
    array = None # no predefined polygon
    min_rand_coord = 1
    max_rand_coord = 10000
    points_num = 30
    
    crt = Create_random_polygon(array, min_rand_coord, max_rand_coord, points_num)
    polygon_array = crt.main(plot = True)    

==========

Answer 3

As someone said, the minimal length solution is exactly the traveling salesman problem. Here's a non-optimal but feasible approach:

Compute a Delauney triangulation of your points. Successively remove boundary segments until you are left with a boundary that interpolates all points or no more segments can be removed. Don't remove boundary segments if all points of the triangle using that segment are on the boundary. Take this boundary as your path.

I implemented this in Mathematica using 40 random points. Here is a typical result:

The obvious objection is that you might get to a point where not all your points are boundary points, but you can't remove a boundary segment without making the boundary self intersecting. This is a valid objection. It took me dozens of runs to see a case where this happened, but finally got this case:

You can probably see some obvious ways of fixing this using the local topology, but I'll leave the details to you! One thing that might help is "edge flipping" where you take two triangles that share a side, say triangle (p,q,r) and (q,p,s) and replace them with (r,p,s) and (r,s,q) (all coordinates counterclockwise around the triangle). This can be done as long as the resulting triangles in this transformation are also counterclockwise ordered.

To reduce the need for fix-ups, you will want to make good choices of the segments to remove at each step. I used the ratio of the length of the boundary segment to the sum of the lengths of the other side of the candidate triangle (the triangle formed by the potentially incoming point with the segment).

Answer 4

I just had this very same problem and came up with some quite simple solution, also of n*log(n) complexity.

First take some point internal to the figure, it does not matter which, it makes sense for it to be the central point, either in the middle of the most distant points or in the average of all points (like a gravity center).

Then sort all the points according to an angle from which they are seen from the central point. The sort key would be equivalent to atan2 for a point and the center.

That's it. Assuming that p is an array of points (x, y), this is the Python code.

center = reduce(lambda a, b: (a[0] + b[0], a[1] + b[1]), p, (0, 0))
center = (center[0] / len(p), (center[1] / len(p)))
p.sort(key = lambda a: math.atan2(a[1] - center[1], a[0] - center[0]))

Answer 5

What you are seeking is called a simple polygonization in the literature. See, for example, this web page on the topic. It is easy to generate a star-shaped polygonization, as Miguel says, but difficult to find, for example, a minimal perimeter polygonization, which is a minimal TSP, as Axel Kemper mentions. There are in general an exponential number of different polygonizations for a given point set.

         

For the star-shaped polygonization, there is one issue that requires some attention: the extra point x (in the "kernel" of the star) must not coincide with an existing point! Here is one algorithm to guarantee x. Find the closest pair of points (a,b), and let x be the midpoint of segment ab. Then proceed as per Miguel's post.

Answer 6

Well, if you don't actually care about minimality or anything like that, you can just place new point inside the convex hull and then order the other points by angle to this new point. You'll get a non-intersecting polygon.

Answer 7

I modified the codes in Comrade Che 's answer to avoid generating intersecting polygon when there exit more than one leftmost or rightmost points(e.g., [(10, 20), (17, 5), (1, 16), (1, 14), (20, 8), (4, 7), (6, 9)]). The main change is that if there exit more than one leftmost or rightmost points, then compare with their y coordinates and select the bottom one as the leftmost or the rightmost point. Here are the codes:

import random
from operator import itemgetter
import numpy
import matplotlib
import matplotlib.pyplot

class Create_random_polygon:

def __init__(self, array, min_rand_coord = None, max_rand_coord = None, points_num = None):        
    self.array = array
    self.min_rand_coord = min_rand_coord 
    self.max_rand_coord = max_rand_coord
    self.points_num = points_num

def generate_random_points(self):
    random_coords_list = []
    for x in range(self.points_num):
        coords_tuple = (random.randint(self.min_rand_coord, self.max_rand_coord),
                        random.randint(self.min_rand_coord, self.max_rand_coord))
        random_coords_list.append(coords_tuple)
    self.array = random_coords_list
    return random_coords_list

def close_line_to_polygon(self):
    a = self.array[0]
    b = self.array[len(self.array)-1]
    if a == b:
        pass
    else:
        self.array.append(a)    

def find_leftmost_point(self):
    leftmost_point = None
    leftmost_x = None
    leftmost_y = None
    for point in self.array:
        x = point[0]
        y = point[1]
        if (leftmost_x == None) or (x < leftmost_x) or (x == leftmost_x and y < leftmost_y):
            leftmost_x = x
            leftmost_y = y
            leftmost_point = point
    return leftmost_point

def find_rightmost_point(self):
    rightmost_point = None
    rightmost_x = None
    rightmost_y = None
    for point in self.array:
        x = point[0]
        y = point[1]
        if (rightmost_x == None) or (x > rightmost_x) or (x == rightmost_x and y < rightmost_y ):
            rightmost_x = x
            rightmost_y = y
            rightmost_point = point
    return rightmost_point

def is_point_above_the_line(self, point, line_points):
    """return 1 if point is above the line
       return -1 if point is below the line
       return  0 if point is lays on the line"""
    px, py = point
    P1, P2 = line_points
    P1x, P1y = P1[0], P1[1]
    P2x, P2y = P2[0], P2[1]
    array = numpy.array([
        [P1x - px, P1y - py],
        [P2x - px, P2y - py],
        ])
    det = numpy.linalg.det(array)
    sign = numpy.sign(det)
    return sign

def sort_array_into_A_B_C(self, line_points):
    [(x_lm, y_lm), (x_rm, y_rm)] = line_points
    A_array, B_array, C_array = [], [], []
    for point in self.array:
        x, y = point
        sing = self.is_point_above_the_line( (x, y), line_points)
        if sing == 0:
            C_array.append(point)
        elif sing == -1:
            A_array.append(point)
        elif sing == 1:
            B_array.append(point)
    return A_array, B_array, C_array

def sort_and_merge_A_B_C_arrays(self, A_array, B_array, C_array):
    A_C_array = [*A_array, *C_array]
    A_C_array.sort(key=itemgetter(0))
    B_array.sort(key=itemgetter(0), reverse=True)        
    merged_arrays = [*A_C_array, *B_array]
    self.array = merged_arrays

def show_image(self, array, line_points, A_array, B_array, C_array):
    [(x_lm, y_lm), (x_rm, y_rm)] = line_points        
    x = [x[0] for x in array]
    y = [y[1] for y in array]
    Ax = [x[0] for x in A_array]
    Ay = [y[1] for y in A_array]
    Bx = [x[0] for x in B_array]
    By = [y[1] for y in B_array]
    Cx = [x[0] for x in C_array]
    Cy = [y[1] for y in C_array]          
    matplotlib.pyplot.plot(Ax, Ay, 'o', c='orange') # below the line
    matplotlib.pyplot.plot(Bx, By, 'o', c='blue') # above the line
    matplotlib.pyplot.plot(Cx, Cy, 'o', c='black') # on the line
    matplotlib.pyplot.plot(x_lm, y_lm, 'o', c='green') # leftmost point
    matplotlib.pyplot.plot(x_rm, y_rm, 'o', c='red') # rightmost point
    x_plot = matplotlib.pyplot.plot([x_lm, x_rm], [y_lm, y_rm], linestyle=':', color='black', linewidth=0.5) # polygon's division line
    x_plot = matplotlib.pyplot.plot(x, y, color='black', linewidth=1) # connect points by line in order of apperiance        
    matplotlib.pyplot.show()

def main(self, plot = False):
    'First output is random polygon coordinates array (other stuff for ploting)'
    print(self.array)
    if self.array == None:
        if not all(
            [isinstance(min_rand_coord, int),
             isinstance(max_rand_coord, int),
             isinstance(points_num, int),]
            ):
            print('Error! Values must be "integer" type:', 'min_rand_coord =',min_rand_coord, ', max_rand_coord =',max_rand_coord, ', points_num =',points_num)
        else:                
            self.array = self.generate_random_points()            

    print(self.array)
    x_lm, y_lm = self.find_leftmost_point()
    x_rm, y_rm = self.find_rightmost_point()
    line_points = [(x_lm, y_lm), (x_rm, y_rm)]

    A_array, B_array, C_array = self.sort_array_into_A_B_C(line_points)
    self.sort_and_merge_A_B_C_arrays(A_array, B_array, C_array)
    self.close_line_to_polygon()
    if plot:
        self.show_image(self.array, line_points, A_array, B_array, C_array)
    return self.array

if __name__ == "__main__":
# predefined polygon
 array = [ 
    (0, 0),
    (2, 2),
    (4, 4),
    (5, 5),
    (0, 5),        
    (1, 4),
    (4, 2),
    (3, 3),
    (2, 1),
    (5, 0),
    ]    
 #array = [(10, 20), (17, 5), (1, 16), (1, 14), (20, 8), (4, 7), (6, 9)]
 #array = [(1, 19), (12, 18), (10, 1), (1, 9), (5, 16), (10, 18), (2, 1)]
 #array = [(13, 17), (15, 3), (14, 13), (11, 8), (7, 16), (7, 7), (10, 15)] 
 array = None # no predefined polygon
 min_rand_coord = 1
 max_rand_coord = 10000
 points_num = 30

 crt = Create_random_polygon(array, min_rand_coord, max_rand_coord, points_num)
 polygon_array = crt.main(plot = True)  

Answer 8

Here is my Typescript implementation of Pawel Pieczul's answer which worked perfectly for my use case involving simple polygons:

interface Point {
    x: number,
    y: number,
    z?: number,
}

const getCentroid = (points: Point[]) => {
    let centroid = { x: 0, y: 0 }
    for (let i = 0; i < points.length; i++) {
        centroid.x += points[i].x
        centroid.y += points[i].y
    }

    centroid.x /= points.length
    centroid.y /= points.length
    return centroid
}

export const sortNonIntersecting = (points: Point[]) => {
    const center = getCentroid(points)
    return points.slice().sort((a: Point, b: Point) => {
        const angleA = Math.atan2(a.y - center.y, a.x - center.x)
        const angleB = Math.atan2(b.y - center.y, b.x - center.x)
        return angleA - angleB
    })
}

Answer 9

I believe you can use Graham scan algorithm to solve your problem.

Edit: in general, Convex hull algorithms are something to look at.

Answer 10

Testing if two segments intersect is simple and fast. See that for example.

With that you could build your polygon iteratively :

Source Points : S = {S0, ... Si, Sj,...}

Final Polygon : A = {A0, ... Ai, Aj,...}

You start with S full and A empty.

Take the first 3 points of S and move them to A. This triangle is of course not self intersecting.

Then, until S is empty, take its first remaining point, that we'll call P, and look for a position in A where it could be inserted.

This position is i+1 for the first i verifying that neither [Ai-P] nor [Ai+1-P] intersects any other segments [Ak-Ak+1].

Your new polygon A is thus {A0, ... Ai, P, Ai+1, ...}

Answer 11

Flutter and Dart developers can use this. I am using this for fix user selected points for create a polygon. When users draw polygons on the map, they generally do not mark points in order.

Example Result: Left one corrected by using this method, Right one is not.

Here is the dart implementation of Pawel's answer;

      LatLng findCentroid(List<LatLng> points) {
        double x = 0;
        double y = 0;
        for (LatLng p in points) {
          x += p.latitude;
          y += p.longitude;
        }
        LatLng center = new LatLng(0, 0);
        center.latitude = x / points.length;
        center.longitude = y / points.length;
        return center;
      }
    
      List<LatLng> sortVerticies(List<LatLng> points) {
        // get centroid
        LatLng center = findCentroid(points);
    
        points.sort((a, b){
          double a1 = (radsToDegrees(Math.atan2(a.latitude - center.latitude, a.longitude - center.longitude)) + 360) % 360;
          double a2 = (radsToDegrees(Math.atan2(b.latitude - center.latitude, b.longitude - center.longitude)) + 360) % 360;
          return (a1 - a2).toInt();
        });
        return points;
      }
    
      num radsToDegrees(num rad) {
        return (rad * 180.0) / Math.pi;
      }

Answer 12

I had a similar question and wasn't able to find a good answer, but finally found this solution that does a great job of explaining it. If you're using Python you can use NumPy for it pretty easily.

import numpy as np


coords = np.array([[12,43], [200,300], [76,170], [260,78], [23,90],[170,323], [54,39], [254, 89]])

def sort_coordinates(list_of_xy_coords):
    cx, cy = list_of_xy_coords.mean(0)
    x, y = list_of_xy_coords.T
    angles = np.arctan2(x-cx, y-cy)
    indices = np.argsort(angles)
    return list_of_xy_coords[indices]