Algorithm to generate random 2D polygon

For a convex 2D polygon (totally off the top of my head):

  1. Generate a random radius, R

  2. Generate N random points on the circumference of a circle of Radius R

  3. Move around the circle and draw straight lines between adjacent points on the circle.


There's a neat way to do what you want by taking advantage of the MATLAB classes DelaunayTri and TriRep and the various methods they employ for handling triangular meshes. The code below follows these steps to create an arbitrary simple polygon:

  • Generate a number of random points equal to the desired number of sides plus a fudge factor. The fudge factor ensures that, regardless of the result of the triangulation, we should have enough facets to be able to trim the triangular mesh down to a polygon with the desired number of sides.

  • Create a Delaunay triangulation of the points, resulting in a convex polygon that is constructed from a series of triangular facets.

  • If the boundary of the triangulation has more edges than desired, pick a random triangular facet on the edge that has a unique vertex (i.e. the triangle only shares one edge with the rest of the triangulation). Removing this triangular facet will reduce the number of boundary edges.

  • If the boundary of the triangulation has fewer edges than desired, or the previous step was unable to find a triangle to remove, pick a random triangular facet on the edge that has only one of its edges on the triangulation boundary. Removing this triangular facet will increase the number of boundary edges.

  • If no triangular facets can be found matching the above criteria, post a warning that a polygon with the desired number of sides couldn't be found and return the x and y coordinates of the current triangulation boundary. Otherwise, keep removing triangular facets until the desired number of edges is met, then return the x and y coordinates of triangulation boundary.

Here's the resulting function:

function [x, y, dt] = simple_polygon(numSides)

    if numSides < 3
        x = [];
        y = [];
        dt = DelaunayTri();
        return
    end

    oldState = warning('off', 'MATLAB:TriRep:PtsNotInTriWarnId');

    fudge = ceil(numSides/10);
    x = rand(numSides+fudge, 1);
    y = rand(numSides+fudge, 1);
    dt = DelaunayTri(x, y);
    boundaryEdges = freeBoundary(dt);
    numEdges = size(boundaryEdges, 1);

    while numEdges ~= numSides
        if numEdges > numSides
            triIndex = vertexAttachments(dt, boundaryEdges(:,1));
            triIndex = triIndex(randperm(numel(triIndex)));
            keep = (cellfun('size', triIndex, 2) ~= 1);
        end
        if (numEdges < numSides) || all(keep)
            triIndex = edgeAttachments(dt, boundaryEdges);
            triIndex = triIndex(randperm(numel(triIndex)));
            triPoints = dt([triIndex{:}], :);
            keep = all(ismember(triPoints, boundaryEdges(:,1)), 2);
        end
        if all(keep)
            warning('Couldn''t achieve desired number of sides!');
            break
        end
        triPoints = dt.Triangulation;
        triPoints(triIndex{find(~keep, 1)}, :) = [];
        dt = TriRep(triPoints, x, y);
        boundaryEdges = freeBoundary(dt);
        numEdges = size(boundaryEdges, 1);
    end

    boundaryEdges = [boundaryEdges(:,1); boundaryEdges(1,1)];
    x = dt.X(boundaryEdges, 1);
    y = dt.X(boundaryEdges, 2);

    warning(oldState);

end

And here are some sample results:

enter image description here

The generated polygons could be either convex or concave, but for larger numbers of desired sides they will almost certainly be concave. The polygons are also generated from points randomly generated within a unit square, so polygons with larger numbers of sides will generally look like they have a "squarish" boundary (such as the lower right example above with the 50-sided polygon). To modify this general bounding shape, you can change the way the initial x and y points are randomly chosen (i.e. from a Gaussian distribution, etc.).


I took @MitchWheat and @templatetypedef's idea of sampling points on a circle and took it a bit farther.

In my application I need to be able to control how weird the polygons are, ie start with regular polygons and as I crank up the parameters they get increasingly chaotic. The basic idea is as stated by @templatetypedef; walk around the circle taking a random angular step each time, and at each step put a point at a random radius. In equations I'm generating the angular steps as equations for the angles and radii of the vertices

where theta_i and r_i give the angle and radius of each point relative to the centre, U(min, max) pulls a random number from a uniform distribution, and N(mu, sigma) pulls a random number from a Gaussian distribution, and clip(x, min, max) thresholds a value into a range. This gives us two really nice parameters to control how wild the polygons are - epsilon which I'll call irregularity controls whether or not the points are uniformly space angularly around the circle, and sigma which I'll call spikeyness which controls how much the points can vary from the circle of radius r_ave. If you set both of these to 0 then you get perfectly regular polygons, if you crank them up then the polygons get crazier.

I whipped this up quickly in python and got stuff like this: some polygons I generated

Here's the full python code:

import math, random
from typing import List, Tuple


def generate_polygon(center: Tuple[float, float], avg_radius: float,
                     irregularity: float, spikiness: float,
                     num_vertices: int) -> List[Tuple[float, float]]:
    """
    Start with the center of the polygon at center, then creates the
    polygon by sampling points on a circle around the center.
    Random noise is added by varying the angular spacing between
    sequential points, and by varying the radial distance of each
    point from the centre.

    Args:
        center (Tuple[float, float]):
            a pair representing the center of the circumference used
            to generate the polygon.
        avg_radius (float):
            the average radius (distance of each generated vertex to
            the center of the circumference) used to generate points
            with a normal distribution.
        irregularity (float):
            variance of the spacing of the angles between consecutive
            vertices.
        spikiness (float):
            variance of the distance of each vertex to the center of
            the circumference.
        num_vertices (int):
            the number of vertices of the polygon.
    Returns:
        List[Tuple[float, float]]: list of vertices, in CCW order.
    """
    # Parameter check
    if irregularity < 0 or irregularity > 1:
        raise ValueError("Irregularity must be between 0 and 1.")
    if spikiness < 0 or spikiness > 1:
        raise ValueError("Spikiness must be between 0 and 1.")

    irregularity *= 2 * math.pi / num_vertices
    spikiness *= avg_radius
    angle_steps = random_angle_steps(num_vertices, irregularity)

    # now generate the points
    points = []
    angle = random.uniform(0, 2 * math.pi)
    for i in range(num_vertices):
        radius = clip(random.gauss(avg_radius, spikiness), 0, 2 * avg_radius)
        point = (center[0] + radius * math.cos(angle),
                 center[1] + radius * math.sin(angle))
        points.append(point)
        angle += angle_steps[i]

    return points
def random_angle_steps(steps: int, irregularity: float) -> List[float]:
    """Generates the division of a circumference in random angles.

    Args:
        steps (int):
            the number of angles to generate.
        irregularity (float):
            variance of the spacing of the angles between consecutive vertices.
    Returns:
        List[float]: the list of the random angles.
    """
    # generate n angle steps
    angles = []
    lower = (2 * math.pi / steps) - irregularity
    upper = (2 * math.pi / steps) + irregularity
    cumsum = 0
    for i in range(steps):
        angle = random.uniform(lower, upper)
        angles.append(angle)
        cumsum += angle

    # normalize the steps so that point 0 and point n+1 are the same
    cumsum /= (2 * math.pi)
    for i in range(steps):
        angles[i] /= cumsum
    return angles
def clip(value, lower, upper):
    """
    Given an interval, values outside the interval are clipped to the interval
    edges.
    """
    return min(upper, max(value, lower))

@MateuszKonieczny here is code to create an image of a polygon from a list of vertices.

vertices = generate_polygon(center=(250, 250),
                            avg_radius=100,
                            irregularity=0.35,
                            spikiness=0.2,
                            num_vertices=16)

black = (0, 0, 0)
white = (255, 255, 255)
img = Image.new('RGB', (500, 500), white)
im_px_access = img.load()
draw = ImageDraw.Draw(img)

# either use .polygon(), if you want to fill the area with a solid colour
draw.polygon(vertices, outline=black, fill=white)

# or .line() if you want to control the line thickness, or use both methods together!
draw.line(vertices + [vertices[0]], width=2, fill=black)

img.show()

# now you can save the image (img), or do whatever else you want with it.

As @templatetypedef and @MitchWheat said, it is easy to do so by generating N random angles and radii. It is important to sort the angles, otherwise it will not be a simple polygon. Note that I am using a neat trick to draw closed curves - I described it in here. By the way, the polygons might be concave.

Note that all of these polygons will be star shaped. Generating a more general polygon is not a simple problem at all. Just to give you a taste of the problem - check out http://www.cosy.sbg.ac.at/~held/projects/rpg/rpg.html and http://compgeom.cs.uiuc.edu/~jeffe/open/randompoly.html.

enter image description here

function CreateRandomPoly()
    figure();
    colors = {'r','g','b','k'};
    for i=1:5
        [x,y]=CreatePoly();
        c = colors{ mod(i-1,numel(colors))+1};
        plotc(x,y,c);
        hold on;
    end        
end

function [x,y]=CreatePoly()
    numOfPoints = randi(30);
    theta = randi(360,[1 numOfPoints]);
    theta = theta * pi / 180;
    theta = sort(theta);
    rho = randi(200,size(theta));
    [x,y] = pol2cart(theta,rho);    
    xCenter = randi([-1000 1000]);
    yCenter = randi([-1000 1000]);
    x = x + xCenter;
    y = y + yCenter;    
end

function plotc(x,y,varargin)
    x = [x(:) ; x(1)];
    y = [y(:) ; y(1)];
    plot(x,y,varargin{:})
end