How can I generate a set of points evenly distribu

2019-02-06 11:39发布

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?

7条回答
做个烂人
2楼-- · 2019-02-06 11:48

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);
            points.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;
    }
查看更多
Fickle 薄情
3楼-- · 2019-02-06 11:53

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()
查看更多
走好不送
4楼-- · 2019-02-06 11:53

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
查看更多
Melony?
5楼-- · 2019-02-06 11:54

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.

查看更多
霸刀☆藐视天下
6楼-- · 2019-02-06 11:56

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;
}

}
查看更多
Summer. ? 凉城
7楼-- · 2019-02-06 12: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.

查看更多
登录 后发表回答