How to generate a random convex polygon?

Question

I'm trying to devise a method for generating random 2D convex polygons. It has to have the following properties:

• coordinates should be integers;
• the polygon should lie inside a square with corners (0, 0) and (C, C), where C is given;
• the polygon should have number of vertices close to a given number N.

For example, generate random polygons that have 10 vertices and lie inside square [0..100]x[0..100].

What makes this task hard, is the fact that the coordinates should be integers.

The approach I tried was to generate random set of points in the given square and compute the convex hull of these points. But the resultant convex hull is very little vertices compared to N.

Any ideas?

Innocency

Answer 1

Here is the fastest algorithm I know that generates each convex polygon with equal probability. The output has exactly N vertices, and the running time is O(N log N), so it can generate even large polygons very quickly.

• Generate two lists, X and Y, of N random integers between 0 and C. Make sure there are no duplicates.
• Sort X and Y and store their maximum and minimum elements.
• Randomly divide the other (not max or min) elements into two groups: X1 and X2, and Y1 and Y2.
• Re-insert the minimum and maximum elements at the start and end of these lists (minX at the start of X1 and X2, maxX at the end, etc.).
• Find the consecutive differences (X1[i + 1] - X1[i]), reversing the order for the second group (X2[i] - X2[i + 1]). Store these in lists XVec and YVec.
• Randomize (shuffle) YVec and treat each pair XVec[i] and YVec[i] as a 2D vector.
• Sort these vectors by angle and then lay them end-to-end to form a polygon.
• Move the polygon to the original min and max coordinates.

An animation and Java implementation is available here: Generating Random Convex Polygons.

This algorithm is based on a paper by Pavel Valtr: “Probability that n random points are in convex position.” Discrete & Computational Geometry 13.1 (1995): 637-643.

Automatism

What are you supposed to do in the case where N == 3 since x1 or x2 will have to be empty then? – Irish

X1 and X2 both include the minimum and maximum element after the next step. The middle element is added to one or the other, randomly. You're essentially deciding whether the middle point is on the top or bottom of the triangle. Of course, all this is overkill for N=3, as any three points form a unique triangle and all triangles are convex. – Automatism

Answer 2

Following @Mangara answer there is JAVA implementation, if someone is interested in Python port of it

import random
from math import atan2


def to_convex_contour(vertices_count,
                      x_generator=random.random,
                      y_generator=random.random):
    """
    Port of Valtr algorithm by Sander Verdonschot.

    Reference:
        http://cglab.ca/~sander/misc/ConvexGeneration/ValtrAlgorithm.java

    >>> contour = to_convex_contour(20)
    >>> len(contour) == 20
    True
    """
    xs = [x_generator() for _ in range(vertices_count)]
    ys = [y_generator() for _ in range(vertices_count)]
    xs = sorted(xs)
    ys = sorted(ys)
    min_x, *xs, max_x = xs
    min_y, *ys, max_y = ys
    vectors_xs = _to_vectors_coordinates(xs, min_x, max_x)
    vectors_ys = _to_vectors_coordinates(ys, min_y, max_y)
    random.shuffle(vectors_ys)

    def to_vector_angle(vector):
        x, y = vector
        return atan2(y, x)

    vectors = sorted(zip(vectors_xs, vectors_ys),
                     key=to_vector_angle)
    point_x = point_y = 0
    min_polygon_x = min_polygon_y = 0
    points = []
    for vector_x, vector_y in vectors:
        points.append((point_x, point_y))
        point_x += vector_x
        point_y += vector_y
        min_polygon_x = min(min_polygon_x, point_x)
        min_polygon_y = min(min_polygon_y, point_y)
    shift_x, shift_y = min_x - min_polygon_x, min_y - min_polygon_y
    return [(point_x + shift_x, point_y + shift_y)
            for point_x, point_y in points]


def _to_vectors_coordinates(coordinates, min_coordinate, max_coordinate):
    last_min = last_max = min_coordinate
    result = []
    for coordinate in coordinates:
        if _to_random_boolean():
            result.append(coordinate - last_min)
            last_min = coordinate
        else:
            result.append(last_max - coordinate)
            last_max = coordinate
    result.extend((max_coordinate - last_min,
                   last_max - max_coordinate))
    return result


def _to_random_boolean():
    return random.getrandbits(1)

Gownsman

Answer 3

This isn't quite complete, but it may give you some ideas.

Bail out if N < 3. Generate a unit circle with N vertices, and rotate it random [0..90] degrees.

Randomly extrude each vertex outward from the origin, and use the sign of the cross product between each pair of adjacent vertices and the origin to determine convexity. This is the step where there are tradeoffs between speed and quality.

After getting your vertices set up, find the vertex with the largest magnitude from the origin. Divide every vertex by that magnitude to normalize the polygon, and then scale it back up by (C/2). Translate to (C/2, C/2) and cast back to integer.

Morale

Here's a link for more info on convexity testing: gamedev.net/topic/… – Morale

This works for floating point coordinates. How do you make sure that when you cast back to integers, the polygon does not become concave? You could remove the problematic vertices, but this could dramatically reduce the number of vertices. – Innocency

Answer 4

A simple algorithm would be:

1. Start with random line (a two vertices and two edges polygon)
2. Take random edge E of the polygon
3. Make new random point P on this edge
4. Take a line L perpendicular to E going through point P. By calculating intersection between line T and lines defined by the two edges adjacent to E, calculate the maximum offset of P when the convexity is not broken.
5. Offset the point P randomly in that range.
6. If not enough points, repeat from 2.

Matriarch

How does that support integer coordinates? – Innocency

@Jasiu: I don't see how could it not support integer coordinates. Just calculate all generated points in integers and maybe clamp the result to the desired range. – Matriarch

Answer 5

I've made the ruby port as well thanks to both @Mangara's answer and @Azat's answer:

#!/usr/bin/env ruby
# frozen_string_literal: true

module ValtrAlgorithm
  module_function def random_polygon(length)
    raise ArgumentError, "length should be > 2" unless length > 2

    min_x, *xs, max_x = Array.new(length) { rand }.sort
    min_y, *ys, max_y = Array.new(length) { rand }.sort
    # Divide the interior points into two chains and
    # extract the vector components.
    vec_xs = to_random_vectors(xs, min_x, max_x)
    vec_ys = to_random_vectors(ys, min_y, max_y).
      # Randomly pair up the X- and Y-components
      shuffle
    # Combine the paired up components into vectors
    vecs = vec_xs.zip(vec_ys).
      # Sort the vectors by angle, in a counter clockwise fashion. Remove the
      # `-` to make it clockwise.
      sort_by { |x, y| - Math.atan2(y, x) }

    # Lay them end-to-end
    point_x = point_y = 0
    min_polygon_x = min_polygon_y = 0
    points = []
    vecs.each do |vec_x, vec_y|
      points.append([vec_x, vec_y])
      point_x += vec_x
      point_y += vec_y
      min_polygon_x = [min_polygon_x, point_x].min
      min_polygon_y = [min_polygon_y, point_y].min
    end
    shift_x = min_x - min_polygon_x
    shift_y = min_y - min_polygon_y
    result = points.map { |point_x, point_y| [point_x + shift_x, point_y + shift_y] }
    # Append first point to make it a valid linear ring
    result << result.first
  end

  private def to_random_vectors(coordinates, min, max)
    last_min = last_max = min
    ary = []
    coordinates.each do |coordinate|
      if rand > 0.5
        ary << coordinate - last_min
        last_min = coordinate
      else
        ary << last_max - coordinate
        last_max = coordinate
      end
    end
    ary << max - last_min << last_max - max
  end
end

Kronos

Answer 6

Here's another version of Valtr's algorithm using numpy. :)

import numpy as np
import numpy.typing and npt


def generateConvex(n: int) -> npt.NDArray[np.float64]:
    '''
    Generate convex shappes according to Pavel Valtr's 1995 alogrithm. Ported from
    Sander Verdonschot's Java version, found here:
    https://cglab.ca/~sander/misc/ConvexGeneration/ValtrAlgorithm.java
    '''
    # initialise random coordinates
    X_rand, Y_rand = np.sort(np.random.random(n)), np.sort(np.random.random(n))
    X_new, Y_new = np.zeros(n), np.zeros(n)

    # divide the interior points into two chains
    last_true = last_false = 0
    for i in range(1, n):
        if i != n - 1:
            if random.getrandbits(1):
                X_new[i] = X_rand[i] - X_rand[last_true]
                Y_new[i] = Y_rand[i] - Y_rand[last_true]
                last_true = i
            else:
                X_new[i] = X_rand[last_false] - X_rand[i]
                Y_new[i] = Y_rand[last_false] - Y_rand[i]
                last_false = i
        else:
            X_new[0] = X_rand[i] - X_rand[last_true]
            Y_new[0] = Y_rand[i] - Y_rand[last_true]
            X_new[i] = X_rand[last_false] - X_rand[i]
            Y_new[i] = Y_rand[last_false] - Y_rand[i]

    # randomly combine x and y and sort by polar angle
    np.random.shuffle(Y_new)
    vertices = np.stack((X_new, Y_new), axis=-1)
    vertices = vertices[np.argsort(np.arctan2(vertices[:, 1], vertices[:, 0]))]

    # arrange points end to end to form a polygon
    vertices = np.cumsum(vertices, axis=0)

    # center around the origin
    x_max, y_max = np.max(vertices[:, 0]), np.max(vertices[:, 1])
    vertices[:, 0] += ((x_max - np.min(vertices[:, 0])) / 2) - x_max
    vertices[:, 1] += ((y_max - np.min(vertices[:, 1])) / 2) - y_max

    return vertices

Jotting

Answer 7

Here is an C++11 realization of the Pavel Valtr algorithm introduced in Mangara's Answer with some tricks similar to lewiswolf's Answer and more randomness realised by separated division process of X and Y coordinate.

#include <algorithm>
#include <iostream>
#include <random>

struct randPoly {
  int RND_MAX = 655369;
  std::random_device dev;
  std::mt19937 rng;
  std::uniform_int_distribution<std::mt19937::result_type> random_numer;
  std::uniform_int_distribution<std::mt19937::result_type> random_logic;

  randPoly() : rng(dev()), random_numer(0, RND_MAX), random_logic(0, 1) {}
  virtual ~randPoly() {}

  int generate(const int n, const double r0, std::vector<double>& poly_x,
               std::vector<double>& poly_y) {
    auto gen = [&]() { return random_numer(rng); };
    // initialize random samples and sort them
    int m = n / 2;
    std::vector<int> x(n), y(n), vx(n), vy(n), idx(n);
    std::vector<double> a(n);
    std::generate(x.begin(), x.end(), gen);
    std::generate(y.begin(), y.end(), gen);
    std::iota(idx.begin(), idx.end(), 0);
    std::sort(x.begin(), x.end());
    std::sort(y.begin(), y.end());
    // divide samples and get vector component
    int x0 = x[0], x1 = x0;
    for (int k = 1; k < n - 1; ++k) {
      if (random_logic(rng)) {
        vx[k - 1] = x[k] - x0;
        x0 = x[k];
      } else {
        vx[k - 1] = x1 - x[k];
        x1 = x[k];
      }
    }
    vx[n - 2] = x[n - 1] - x0;
    vx[n - 1] = x1 - x[n - 1];
    int y0 = y[0], y1 = y0;
    for (int k = 1; k < n - 1; ++k) {
      if (random_logic(rng)) {
        vy[k - 1] = y[k] - y0;
        y0 = y[k];
      } else {
        vy[k - 1] = y1 - y[k];
        y1 = y[k];
      }
    }
    vy[n - 2] = y[n - 1] - y0;
    vy[n - 1] = y1 - y[n - 1];
    // random pair up vector components and sort by angle
    std::shuffle(vy.begin(), vy.end(), rng);
    for (int k = 0; k < n; ++k) {
      a[k] = std::atan2(vy[k], vx[k]);
    }
    std::sort(idx.begin(), idx.end(),
              [&a](int& lhs, int& rhs) { return a[lhs] < a[rhs]; });
    // form the polygon by connencting vectors
    double x_max = 0, y_max = 0, x_min = 0, y_min = 0;
    x[0] = y[0] = 0;
    for (int k = 1; k < n; ++k) {
      x[k] = x[k - 1] + vx[idx[k - 1]];
      y[k] = y[k - 1] + vy[idx[k - 1]];
      if (x[k] > x_max) {
        x_max = x[k];
      } else if (x[k] < x_min) {
        x_min = x[k];
      }
      if (y[k] > y_max) {
        y_max = y[k];
      } else if (y[k] < y_min) {
        y_min = y[k];
      }
    }
    // center and resize the polygon
    poly_x.resize(n);
    poly_y.resize(n);
    double x_offset = -(x_max + x_min) / 2.0;
    double y_offset = -(y_max + y_min) / 2.0;
    double scale = r0 / std::max(x_max - x_min, y_max - y_min);
    for (int k = 0; k < n; ++k) {
      poly_x[k] = scale * (x[k] + x_offset);
      poly_y[k] = scale * (y[k] + y_offset);
    }
    return 0;
  }
};

int main(int, char**) {
  randPoly rp;
  std::vector<double> poly_x, poly_y;
  rp.generate(8, 2.0, poly_x, poly_y);
  for (int k = 0; k < poly_x.size(); ++k) {
    std::cout << poly_x[k] << "  " << poly_y[k] << std::endl;
  }
}

Example shown in Rviz

Draught

Answer 8

Your initial approach is correct - calculating the convex hull is the only way you will satisfy randomness, convexity and integerness.

The only way I can think of optimizing your algorithm to get "more points" out is by organizing them around a circle instead of completely randomly. Your points should more likely be near the "edges" of your square than near the center. At the center, the probability should be ~0, since the polygon must be convex.

One simple option would be setting a minimum radius for your points to appear - maybe C/2 or C*0.75. Calculate the center of the C square, and if a point is too close, move it away from the center until a minimum distance is reached.

Husha

Answer 9

I think the above suggestion with Numpy is excellent but the generateConvex() function may be improved by avoiding the for loop. And I would like to add a few lines of code with Matplotlib to visualize the result.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
rng = np.random.default_rng()

def generateConvex(n):    
    # initialise random coordinates
    XY_rand = np.sort(rng.random((n, 2)), axis=0)
    
    # divide the interior points into two chains
    rand_bool = rng.choice([True, False], n-2)
    pos, neg = XY_rand[1:-1][rand_bool], XY_rand[1:-1][~rand_bool]
    
    pos = np.vstack((XY_rand[0], pos, XY_rand[-1]))
    neg = np.vstack((XY_rand[0], neg, XY_rand[-1]))
    vertices = np.vstack((pos[1:] - pos[:-1], neg[:-1] - neg[1:]))
     
    # randomly combine x and y and sort by polar angle
    rng.shuffle(vertices[:,1])
    vertices = vertices[np.argsort(np.arctan2(vertices[:, 1], vertices[:, 0]))]
    
    # arrange points end to end to form a polygon
    vertices = np.cumsum(vertices, axis=0)
    
    # center around the origin
    x_max, y_max = np.max(vertices[:, 0]), np.max(vertices[:, 1])
    vertices[:, 0] -= (x_max + np.min(vertices[:, 0])) / 2
    vertices[:, 1] -= (y_max + np.min(vertices[:, 1])) / 2
    
    return vertices

if __name__ == '__main__':
    n = 42
    p = Polygon(generateConvex(n))
    fig, ax = plt.subplots()

    ax.add_patch(p)
    plt.title(f'{n}-sided convex polygon')
    plt.axis('equal')
    plt.show()

Radley