Uniformly distribute x points inside a circle

Question

I would like to uniformly distribute a predetermined set of points within a circle. By uniform distribution, I mean they should all be equally distanced from each other (hence a random approach won't work). I tried a hexagonal approach, but I had problems consistently reaching the outermost radius.

My current approach is a nested for loop where each outer iteration reduces the radius & number of points, and each inner loop evenly drops points on the new radius. Essentially, it's a bunch of nested circles. Unfortunately, it's far from even. Any tips on how to do this correctly?

Muir

Take a look at low-discrepenacy sequences. – Erotica

What do you want to happen on the boundary, do the points have to be uniformly spaced around the boundary of the circle? Otherwise you could take a uniform grid (triangular, hexagonal, or square) and keep from it only the points within the circle. – Hilariohilarious

RobertDodier, thanks, unfortunately even subrandom numbers don't produce favorable results as the chance of 2 points landing closely next to each other is still relatively high. FamousBlueRaincoat, ideally, yes for this case I would like the outermost points to be uniformly spaced around the edge. – Muir

Sunflower seed arrangement uses golden ratio for angle increments and square root for radius increment. The result looks pretty good, but the outer edge isn't perfectly circular. – Hilariohilarious

Equally spaced points form an hexagonal lattice and can't lie on a circular boundary (except N=7). What you are asking is impossible, you must relax some condition. – Songwriter

Answer 1

The goals of having a uniform distribution within the area and a uniform distribution on the boundary conflict; any solution will be a compromise between the two. I augmented the sunflower seed arrangement with an additional parameter alpha that indicates how much one cares about the evenness of boundary.

alpha=0 gives the typical sunflower arrangement, with jagged boundary:

With alpha=2 the boundary is smoother:

(Increasing alpha further is problematic: Too many points end up on the boundary).

The algorithm places n points, of which the kth point is put at distance sqrt(k-1/2) from the boundary (index begins with k=1), and with polar angle 2*pi*k/phi^2 where phi is the golden ratio. Exception: the last alpha*sqrt(n) points are placed on the outer boundary of the circle, and the polar radius of other points is scaled to account for that. This computation of the polar radius is done in the function radius.

It is coded in MATLAB.

function sunflower(n, alpha)   %  example: n=500, alpha=2
    clf
    hold on
    b = round(alpha*sqrt(n));      % number of boundary points
    phi = (sqrt(5)+1)/2;           % golden ratio
    for k=1:n
        r = radius(k,n,b);
        theta = 2*pi*k/phi^2;
        plot(r*cos(theta), r*sin(theta), 'r*');
    end
end

function r = radius(k,n,b)
    if k>n-b
        r = 1;            % put on the boundary
    else
        r = sqrt(k-1/2)/sqrt(n-(b+1)/2);     % apply square root
    end
end

Hilariohilarious

This works great. For future readers,I made this geodesic by turning theta into a bearing so I could pass it into a destination function. I did that with the following line theta = k * (720 - 360 * phi). It looks correct, but if I'm in error please correct me. – Muir

I think 720 is redundant; when multiplied by k, it's just 2*k full turns. You could use theta = -k*360*phi instead; or theta = k*360*phi since the minus sign is just mirror reflection. – Hilariohilarious

Right again! Looks like I need a refresher in geometry :-/. – Muir

Any idea how to take this into 3D? Plotting all points within a sphere? – Christophe

Answer 2

Might as well tag on my Python translation.

from math import sqrt, sin, cos, pi
phi = (1 + sqrt(5)) / 2  # golden ratio

def sunflower(n, alpha=0, geodesic=False):
    points = []
    angle_stride = 360 * phi if geodesic else 2 * pi / phi ** 2
    b = round(alpha * sqrt(n))  # number of boundary points
    for k in range(1, n + 1):
        r = radius(k, n, b)
        theta = k * angle_stride
        points.append((r * cos(theta), r * sin(theta)))
    return points

def radius(k, n, b):
    if k > n - b:
        return 1.0
    else:
        return sqrt(k - 0.5) / sqrt(n - (b + 1) / 2)


# example
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    fig, ax = plt.subplots()
    points = sunflower(500, alpha=2, geodesic=False)
    xs = [point[0] for point in points]
    ys = [point[1] for point in points]
    ax.scatter(xs, ys)
    ax.set_aspect('equal') # display as square plot with equal axes
    plt.show()

Asquith

Answer 3

Building on top of @OlivelsAWord , here is a Python implementation using numpy:

import numpy as np
import matplotlib.pyplot as plt


def sunflower(n: int, alpha: float) -> np.ndarray:
    # Number of points respectively on the boundary and inside the cirlce.
    n_exterior = np.round(alpha * np.sqrt(n)).astype(int)
    n_interior = n - n_exterior

    # Ensure there are still some points in the inside...
    if n_interior < 1:
        raise RuntimeError(f"Parameter 'alpha' is too large ({alpha}), all "
                           f"points would end-up on the boundary.")
    # Generate the angles. The factor k_theta corresponds to 2*pi/phi^2.
    k_theta = np.pi * (3 - np.sqrt(5))
    angles = np.linspace(k_theta, k_theta * n, n)

    # Generate the radii.
    r_interior = np.sqrt(np.linspace(0, 1, n_interior))
    r_exterior = np.ones((n_exterior,))
    r = np.concatenate((r_interior, r_exterior))

    # Return Cartesian coordinates from polar ones.
    return r * np.stack((np.cos(angles), np.sin(angles)))

    # NOTE: say the returned array is called s. The layout is such that s[0,:]
    # contains X values and s[1,:] contains Y values. Change the above to
    #   return r.reshape(n, 1) * np.stack((np.cos(angles), np.sin(angles)), axis=1)
    # if you want s[:,0] and s[:,1] to contain X and Y values instead.


if __name__ == '__main__':
    fig, ax = plt.subplots()

    # Let's plot three sunflowers with different values of alpha!
    for alpha in (0, 1, 2):
        s = sunflower(500, alpha)
        # NOTE: the 'alpha=0.5' parameter is to control transparency, it does
        # not have anything to do with the alpha used in 'sunflower' ;)
        ax.scatter(s[0], s[1], alpha=0.5, label=f"alpha={alpha}")
    
    # Display as square plot with equal axes and add a legend. Then show the result :)
    ax.set_aspect('equal')
    ax.legend()
    plt.show()

Gluttony

Answer 4

Stumbled across this question and the answer above (so all cred to user3717023 & Matt).
Just adding my translation into R here, in case someone else needed that :)

library(tibble)
library(dplyr)
library(ggplot2)

sunflower <- function(n, alpha = 2, geometry = c('planar','geodesic')) {
  b <- round(alpha*sqrt(n))  # number of boundary points
  phi <- (sqrt(5)+1)/2  # golden ratio

  r <- radius(1:n,n,b)
  theta <- 1:n * ifelse(geometry[1] == 'geodesic', 360*phi, 2*pi/phi^2)

  tibble(
    x = r*cos(theta),
    y = r*sin(theta)
  )
}

radius <- function(k,n,b) {
  ifelse(
    k > n-b,
    1,
    sqrt(k-1/2)/sqrt(n-(b+1)/2)
  )
}

# example:
sunflower(500, 2, 'planar') %>%
    ggplot(aes(x,y)) +
    geom_point()

Bacon

Answer 5

Adding my Java implementation of previous answers with an example (Processing).

int n = 2000; // count of nodes
Float alpha = 2.; // constant that can be adjusted to vary the geometry of points at the boundary
ArrayList<PVector> vertices = new ArrayList<PVector>();
Float scaleFactor = 200.; // scale points beyond their 0.0-1.0 range for visualisation;

void setup() {
  size(500, 500);
  // Test
  vertices = sunflower(n, alpha);
  displayTest(vertices, scaleFactor);
}


ArrayList<PVector> sunflower(int n, Float alpha) {
  Double phi = (1 + Math.sqrt(5)) / 2; // golden ratio
  Double angle = 2 * PI / Math.pow(phi, 2); // value used to calculate theta for each point
  ArrayList<PVector> points = new ArrayList<PVector>();
  Long b = Math.round(alpha*Math.sqrt(n)); // number of boundary points
  Float theta, r, x, y;
  
  for (int i = 1; i < n + 1; i++) {
    r = radius(i, n, b.floatValue());
    theta = i * angle.floatValue();
    x = r * cos(theta);
    y = r * sin(theta);
    PVector p = new PVector(x, y);
    points.add(p);
  }
  return points;
}


Float radius(int k, int n, Float b) {
  if (k > n - b) {
    return 1.0;
  } else {
    Double r = Math.sqrt(k - 0.5) / Math.sqrt(n - (b+1) / 2);
    return r.floatValue();
  }
}


void displayTest(ArrayList<PVector> points, Float size) {
  for (int i = 0; i < points.size(); i++) {
    Float x = size * points.get(i).x;
    Float y = size * points.get(i).y;
    pushMatrix();
    translate(width / 2, height / 2);
    ellipse(x, y, 5, 5);
    popMatrix();
  }
}

Greenbrier

Answer 6

Here's my Unity implementation.

Vector2[] Sunflower(int n, float alpha = 0, bool geodesic = false){
    float phi = (1 + Mathf.Sqrt(5)) / 2;//golden ratio
    float angle_stride = 360 * phi;
    float radius(float k, float n, float b)
    {
        return k > n - b ? 1 : Mathf.Sqrt(k - 0.5f) / Mathf.Sqrt(n - (b + 1) / 2);
    }
    
    int b = (int)(alpha * Mathf.Sqrt(n));  //# number of boundary points

    List<Vector2>points = new List<Vector2>();
    for (int k = 0; k < n; k++)
    {
        float r = radius(k, n, b);
        float theta = geodesic ? k * 360 * phi : k * angle_stride;
        float x = !float.IsNaN(r * Mathf.Cos(theta)) ? r * Mathf.Cos(theta) : 0;
        float y = !float.IsNaN(r * Mathf.Sin(theta)) ? r * Mathf.Sin(theta) : 0;
        points.Add(new Vector2(x, y));
    }   
    return points.ToArray();
}

Cardie