如何生成沿椭圆周长均匀分布的一组点？

Question

如果我想生成一堆均匀分布在一个圆圈周围的点，我可以这样做（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))
]

但是，相同的逻辑不会在椭圆上生成均匀的点：“端”上的点比“边”上的点间距更近。

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

有没有一种简单的方法可以在椭圆周围生成等距点？

Answer 1

这是一个旧线程，但由于我正在寻找沿和椭圆创建均匀间隔点的相同任务并且无法找到实现，因此我提供了实现霍华德伪代码的 Java 代码：

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

}

Answer 2

（已更新：反映新包装）。

可以在FlyingCircus Python 包的数字分支FlyingCircus-Numeric 中找到适用于Python 的有效解决方案。

免责声明：我是它们的主要作者。

简而言之，（简化的）代码看起来（其中a 是短轴，b 是长轴）：

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

它利用scipy.special.ellipeinc() 提供沿椭圆周边的数值积分，以及scipy.optimize.root() 用于求解角度的等弧长度方程。

测试它是否真的有效：

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

Answer 3

您必须计算周长，然后将其分成等长的弧。椭圆的弧长是椭圆积分，不能写成封闭形式，所以需要数值计算。

wolfram 上ellipses 上的文章为您提供了执行此操作所需的公式，但这会很难看。

Answer 4

可能的（数值）计算如下所示：

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

这是一个简单的数值积分方案。如果您需要更高的准确性，您还可以使用任何其他集成方法。

Answer 5

我确定这个帖子早就死了，但我刚刚遇到这个问题，这是最接近解决方案的。

我从 Dave 的回答开始，但我注意到它并没有真正回答发帖人的问题。它不是用弧长等分椭圆，而是用角度。

无论如何，我对他的（很棒的）工作进行了一些调整，以使椭圆除以弧长（这次是用 C# 编写的）。如果你看一下代码，你会看到一些相同的东西 -

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

Answer 6

来自我在 BSE here 的回答。

我将它添加到 stackoverflow 中，因为它是一种不同的方法，它不依赖于固定的迭代步骤，而是依赖于点之间的距离收敛到平均距离。

因此计算时间更短，因为它仅取决于所需的顶点数量和要达到的精度（大约 6 次迭代，不到 0.01%）。

原理是：

0/ 第一步：通常使用 a * cos(t) 和 b * sin(t) 计算点

1/ 计算顶点之间的长度

2/ 根据每个距离与平均距离之间的差距调整角度变化

3/ 重新定位点

4/ 达到所需精度时退出或返回 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()

Answer 7

如果椭圆被压扁，请考虑椭圆周长的公式。 （如果短轴是长轴的三倍）

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

Ellipse Perimeter

Answer 8

有工作的 MATLAB 代码 available here。我在下面复制它，以防该链接失效。学分归原作者所有。

此代码假定主轴是从(x1, y1) 到(x2, y2) 的线段，e 是椭圆的eccentricity。

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