How can I generate a set of points evenly distributed along the perimeter of an ellipse?

Question

If I want to generate a bunch of points distributed uniformly around a circle, I can do this (python):

r = 5  #radius
n = 20 #points to generate
circlePoints = [
    (r * math.cos(theta), r * math.sin(theta))
    for theta in (math.pi*2 * i/n for i in range(n))
]

However, the same logic doesn't generate uniform points on an ellipse: points on the "ends" are more closely spaced than points on the "sides".

r1 = 5
r2 = 10
n = 20 #points to generate
ellipsePoints = [
    (r1 * math.cos(theta), r2 * math.sin(theta))
    for theta in (math.pi*2 * i/n for i in range(n))
]

Is there an easy way to generate equally spaced points around an ellipse?

Answer 1

(UPDATED: to reflect new packaging).

An efficient solution of this problem for Python can be found in the numeric branch FlyingCircus-Numeric, derivated from the FlyingCircus Python package.

Disclaimer: I am the main author of them.

Briefly, the (simplified) code looks (where a is the minor axis, and b is the major axis):

import numpy as np
import scipy as sp
import scipy.optimize

def angles_in_ellipse(
        num,
        a,
        b):
    assert(num > 0)
    assert(a < b)
    angles = 2 * np.pi * np.arange(num) / num
    if a != b:
        e2 = (1.0 - a ** 2.0 / b ** 2.0)
        tot_size = sp.special.ellipeinc(2.0 * np.pi, e2)
        arc_size = tot_size / num
        arcs = np.arange(num) * arc_size
        res = sp.optimize.root(
            lambda x: (sp.special.ellipeinc(x, e2) - arcs), angles)
        angles = res.x 
    return angles

It makes use of scipy.special.ellipeinc() which provides the numerical integral along the perimeter of the ellipse, and scipy.optimize.root() for solving the equal-arcs length equation for the angles.

To test that it is actually working:

a = 10
b = 20
n = 16

phi = angles_in_ellipse(n, a, b)
print(np.round(np.rad2deg(phi), 2))
# [  0.    17.55  36.47  59.13  90.   120.87 143.53 162.45 180.   197.55
#  216.47 239.13 270.   300.87 323.53 342.45]

e = (1.0 - a ** 2.0 / b ** 2.0) ** 0.5
arcs = sp.special.ellipeinc(phi, e)
print(np.round(np.diff(arcs), 4))
# [0.3022 0.2982 0.2855 0.2455 0.2455 0.2855 0.2982 0.3022 0.3022 0.2982
#  0.2855 0.2455 0.2455 0.2855 0.2982]

# plotting
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.gca()
ax.axes.set_aspect('equal')
ax.scatter(b * np.sin(phi), a * np.cos(phi))
plt.show()