20

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?

2
  • 1
    Do you need an accurate description? In this case you'll have to invert the elliptic integral of second kind.
    – Howard
    Aug 7, 2011 at 11:07
  • 1
    C++ has elliptic integral functions starting from C++17 or also some technical report (tr).
    – Pablo H
    Apr 8, 2020 at 21:34

8 Answers 8

17

This is an old thread, but since I am seeking the same task of creating evenly spaced points along and ellipse and was not able to find an implementation, I offer this Java code that implements the pseudo code of Howard:

 package com.math;

  public class CalculatePoints {

  public static void main(String[] args) {
    // TODO Auto-generated method stub

    /*
     * 
        dp(t) = sqrt( (r1*sin(t))^2 + (r2*cos(t))^2)
        circ = sum(dp(t), t=0..2*Pi step 0.0001)

        n = 20

        nextPoint = 0
        run = 0.0
        for t=0..2*Pi step 0.0001
            if n*run/circ >= nextPoint then
                set point (r1*cos(t), r2*sin(t))
                nextPoint = nextPoint + 1
            next
            run = run + dp(t)
        next
     */


    double r1 = 20.0;
    double r2 = 10.0;

    double theta = 0.0;
    double twoPi = Math.PI*2.0;
    double deltaTheta = 0.0001;
    double numIntegrals = Math.round(twoPi/deltaTheta);
    double circ=0.0;
    double dpt=0.0;

    /* integrate over the elipse to get the circumference */
    for( int i=0; i < numIntegrals; i++ ) {
        theta += i*deltaTheta;
        dpt = computeDpt( r1, r2, theta);
        circ += dpt;
    }
    System.out.println( "circumference = " + circ );

    int n=20;
    int nextPoint = 0;
    double run = 0.0;
    theta = 0.0;

    for( int i=0; i < numIntegrals; i++ ) {
        theta += deltaTheta;
        double subIntegral = n*run/circ;
        if( (int) subIntegral >= nextPoint ) {
            double x = r1 * Math.cos(theta);
            double y = r2 * Math.sin(theta);
            System.out.println( "x=" + Math.round(x) + ", y=" + Math.round(y));
            nextPoint++;
        }
        run += computeDpt(r1, r2, theta);
    }
}

static double computeDpt( double r1, double r2, double theta ) {
    double dp=0.0;

    double dpt_sin = Math.pow(r1*Math.sin(theta), 2.0);
    double dpt_cos = Math.pow( r2*Math.cos(theta), 2.0);
    dp = Math.sqrt(dpt_sin + dpt_cos);

    return dp;
}

}
4
  • This worked flawlessly for me. Takes about 160ms on a Droid Maxx (XT1080) running Android 4.4.4. The ellipse uses the full 1280/720 screen with 26 points (one for each letter of the alphabet). I have searched far and wide for a solution and this one is it. Feb 6, 2015 at 3:34
  • This works well for me too. Can you explain what the variables run and dpt mean? Jan 27, 2018 at 2:14
  • 1
    @JeffThompson I'm a few years late here, but I had to squint at it for a while since nobody ever answered your question. "dpt" means dp(t), which calculates the distance from the center of the ellipse to the point on the ellipse at the given theta plus π/2. It's used as an integral term here and summed for a total "size" of the ellipse (not a circumference term at all). "run" is just the current "size" traversed thus far on the second pass over the ellipse (again, offset by π/2). Whenever it passes a multiple of (circ / numpoints), it's time to add a point.
    – Rezer
    Jan 18 at 4:03
  • 1
    The π/2 offset is what makes the points nearly equidistant, it's allocating more points when near the minor axis and less points when near the major axis. To allocate points in a manner that reduces error, it's probably best to just switch the sine and cosine terms in computeDpt which will eliminate the offset and allocate points to where the curve is changing the fastest.
    – Rezer
    Jan 18 at 4:05
11

(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()

plot

10

You have to calculate the perimeter, then divide it into equal length arcs. The length of an arc of an ellipse is an elliptic integral and cannot be written in closed form so you need numerical computation.

The article on ellipses on wolfram gives you the formula needed to do this, but this is going to be ugly.

1
  • 3
    More precisely, for an ellipse (x/a)^2 + (y/b)^2 = 1, the arc length is given by b E(t, e) where E is the incomplete elliptic function of the second kind, e is the eccentricity, and where the ellipse is parametered by (a cos t, b sin t). Aug 9, 2011 at 15:11
5

A possible (numerical) calculation can look as follows:

dp(t) = sqrt( (r1*sin(t))^2 + (r2*cos(t))^2)
circ = sum(dp(t), t=0..2*Pi step 0.0001)

n = 20

nextPoint = 0
run = 0.0
for t=0..2*Pi step 0.0001
    if n*run/circ >= nextPoint then
        set point (r1*cos(t), r2*sin(t))
        nextPoint = nextPoint + 1
    next
    run = run + dp(t)
next

This is a simple numerical integration scheme. If you need better accuracy you might also use any other integration method.

5

I'm sure this thread is long dead by now, but I just came across this issue and this was the closest that came to a solution.

I started with Dave's answer here, but I noticed that it wasn't really answering the poster's question. It wasn't dividing the ellipse equally by arc lengths, but by angle.

Anyway, I made some adjustments to his (awesome) work to get the ellipse to divide equally by arc length instead (written in C# this time). If you look at the code, you'll see some of the same stuff -

    void main()
    {
        List<Point> pointsInEllipse = new List<Point>();
    
        // Distance in radians between angles measured on the ellipse
        double deltaAngle = 0.001;
        double circumference = GetLengthOfEllipse(deltaAngle);
        
        double arcLength = 0.1;
            
        double angle = 0;
            
        // Loop until we get all the points out of the ellipse
        for (int numPoints = 0; numPoints < circumference / arcLength; numPoints++)
        {
            angle = GetAngleForArcLengthRecursively(0, arcLength, angle, deltaAngle);
                
            double x = r1 * Math.Cos(angle);
            double y = r2 * Math.Sin(angle);
            pointsInEllipse.Add(new Point(x, y));
        }
    }

    private double GetLengthOfEllipse()
    {
        // Distance in radians between angles
        double deltaAngle = 0.001;
        double numIntegrals = Math.Round(Math.PI * 2.0 / deltaAngle);

        double radiusX = (rectangleRight - rectangleLeft) / 2;
        double radiusY = (rectangleBottom - rectangleTop) / 2;

        // integrate over the elipse to get the circumference
        for (int i = 0; i < numIntegrals; i++)
        {
            length += ComputeArcOverAngle(radiusX, radiusY, i * deltaAngle, deltaAngle);
        }

        return length;
    }

    private double GetAngleForArcLengthRecursively(double currentArcPos, double goalArcPos, double angle, double angleSeg)
    {

        // Calculate arc length at new angle
        double nextSegLength = ComputeArcOverAngle(majorRadius, minorRadius, angle + angleSeg, angleSeg);

        // If we've overshot, reduce the delta angle and try again
        if (currentArcPos + nextSegLength > goalArcPos) {
            return GetAngleForArcLengthRecursively(currentArcPos, goalArcPos, angle, angleSeg / 2);

            // We're below the our goal value but not in range (
        } else if (currentArcPos + nextSegLength < goalArcPos - ((goalArcPos - currentArcPos) * ARC_ACCURACY)) {
            return GetAngleForArcLengthRecursively(currentArcPos + nextSegLength, goalArcPos, angle + angleSeg, angleSeg);

            // current arc length is in range (within error), so return the angle
        } else
            return angle;
    }

    private double ComputeArcOverAngle(double r1, double r2, double angle, double angleSeg)
    {
        double distance = 0.0;

        double dpt_sin = Math.Pow(r1 * Math.Sin(angle), 2.0);
        double dpt_cos = Math.Pow(r2 * Math.Cos(angle), 2.0);
        distance = Math.Sqrt(dpt_sin + dpt_cos);

        // Scale the value of distance
        return distance * angleSeg;
    }
1
0

From my answer in BSE here .

I add it in stackoverflow as it is a different approach which does not rely on a fixed iteration steps but rely on a convergence of the distances between the points, to the mean distance.

So the calculation is shorter as it depends only on the wanted vertices amount and on the precision to reach (about 6 iterations for less than 0.01%).

The principle is :

0/ First step : calculate the points normally using a * cos(t) and b * sin(t)

1/ Calculate the lengths between vertices

2/ Adjust the angles variations depending on the gap between each distance to the mean distance

3/ Reposition the points

4/ Exit when the wanted precision is reached or return to 1/

import bpy, bmesh
from math import radians, sqrt, cos, sin

rad90 = radians( 90.0 )
rad180 = radians( 180.0 )

def createVertex( bm, x, y ): #uses bmesh to create a vertex
    return bm.verts.new( [x, y, 0] )

def listSum( list, index ): #helper to sum on a list
    sum = 0
    for i in list:
        sum = sum + i[index]
    return sum

def calcLength( points ): #calculate the lenghts for consecutives points
    prevPoint = points[0]
    for point in points :
        dx = point[0] - prevPoint[0]
        dy = point[1] - prevPoint[1]
        dist = sqrt( dx * dx + dy *dy )
        point[3] = dist
        prevPoint = point

def calcPos( points, a, b ): #calculate the positions following the angles
    angle = 0
    for i in range( 1, len(points) - 1 ):
        point = points[i]
        angle += point[2]
        point[0] = a * cos( angle )
        point[1] = b * sin( angle )

def adjust( points ): #adjust the angle by comparing each length to the mean length
    totalLength = listSum( points, 3 )
    averageLength = totalLength / (len(points) - 1)

    maxRatio = 0
    for i in range( 1, len(points) ):
        point = points[i]
        ratio = (averageLength - point[3]) / averageLength
        point[2] = (1.0 + ratio) * point[2]
        absRatio = abs( ratio )
        if absRatio > maxRatio:
            maxRatio = absRatio
    return maxRatio

def ellipse( bm, a, b, steps, limit ):

    delta = rad90 / steps

    angle = 0.0

    points = [] #will be a list of [ [x, y, angle, length], ...]
    for step in range( steps  + 1 ) :
        x = a * cos( angle )
        y = b * sin( angle )
        points.append( [x, y, delta, 0.0] )
        angle += delta

    print( 'start' )
    doContinue = True
    while doContinue:
        calcLength( points )
        maxRatio = adjust( points )
        calcPos( points, a, b )

        doContinue = maxRatio > limit
        print( maxRatio )

    verts = []    
    for point in points:
        verts.append( createVertex( bm, point[0], point[1] ) )

    for i in range( 1, len(verts) ):
        bm.edges.new( [verts[i - 1], verts[i]] )



A = 4
B = 6

bm = bmesh.new()

ellipse( bm, A, B, 32, 0.00001 )

mesh = bpy.context.object.data      
bm.to_mesh(mesh)
mesh.update()
0

Do take into consideration the formula for ellipse perimeter as under if the ellipse is squashed. (If the minor axis is three times as small as the major axis)

tot_size = np.pi*(3*(a+b) -np.sqrt((3*a+b)*a+3*b))

Ellipse Perimeter

-2

There is working MATLAB code available here. I replicate that below in case that link ever goes dead. Credits are due to the original author.

This code assumes that the major axis is a line segment from (x1, y1) to (x2, y2) and e is the eccentricity of the ellipse.

a = 1/2*sqrt((x2-x1)^2+(y2-y1)^2);
b = a*sqrt(1-e^2);
t = linspace(0,2*pi, 20);
X = a*cos(t);
Y = b*sin(t);
w = atan2(y2-y1,x2-x1);
x = (x1+x2)/2 + X*cos(w) - Y*sin(w);
y = (y1+y2)/2 + X*sin(w) + Y*cos(w);
plot(x,y,'o')
axis equal
3
  • This fails the "evenly distributed" requirement, which was the core of the question
    – Eric
    Nov 2, 2016 at 23:22
  • @Eric I'm not sure how this fails that requirement... t is linearly sampled from 0 to 2*pi... Can you elaborate?
    – Maghoumi
    Nov 2, 2016 at 23:47
  • 1
    t might be evenly sampled, but neither the distance nor the angle between points is. This code is equivalent to my example of what I did not want
    – Eric
    Nov 3, 2016 at 13:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.