How to generate a random convex polygon?

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.

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)

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

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

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