Как я могу создать набор точек, равномерно распределенных по периметру эллипса?

Ответ 1

Эффективное решение этой проблемы для Python можно найти в FlyingCircus 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:
        e = (1.0 - a ** 2.0 / b ** 2.0) ** 0.5
        tot_size = sp.special.ellipeinc(2.0 * np.pi, e)
        arc_size = tot_size / num
        arcs = np.arange(num) * arc_size
        res = sp.optimize.root(
            lambda x: (sp.special.ellipeinc(x, e) - 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.    16.4   34.12  55.68  90.   124.32 145.88 163.6  180.   196.4 214.12 235.68 270.   304.32 325.88 343.6 ]

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.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829 0.2829]

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

Ответ 2

Это старый поток, но поскольку я ищу ту же задачу создания равномерно расположенных точек и эллипса и не смог найти реализацию, я предлагаю этот 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;
}

}

Ответ 3

Вы должны вычислить периметр, затем разделите его на равные длины. Длина дуги эллипса является эллиптическим интегралом и не может быть записана в замкнутом виде, поэтому вам нужны численные вычисления.

Статья о эллипсах на вольфраме дает вам формулу, необходимую для этого, но это будет уродливо.

Ответ 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

Это простая схема численного интегрирования. Если вам нужна более высокая точность, вы также можете использовать любой другой метод интеграции.

Ответ 5

Я уверен, что эта ветка уже давно мертва, но я просто столкнулся с этой проблемой, и это было самым близким к решению.

Я начал с ответа Дейва здесь, но я заметил, что это действительно не отвечает на вопрос о плакате. Он не делят эллипс одинаково на длину дуги, а на угол.

Во всяком случае, я внес некоторые корректировки в его (удивительную) работу, чтобы получить эллипс, чтобы разделить его по дуге (вместо этого написано на С#). Если вы посмотрите на код, вы увидите некоторые из тех же материалов -

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

Ответ 6

Из моего ответа в BSE здесь.

Я добавляю его в stackoverflow, поскольку это другой подход, который не полагается на фиксированные шаги итерации, но полагается на конвергенцию расстояний между точками до среднего расстояния.

Таким образом, расчет короче, поскольку он зависит только от количества желаемых вершин и от точности достижения (около 6 итераций менее 0,01%).

Принцип:

0/Первый шаг: вычислить точки, обычно используя * 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()

Ответ 7

Примите во внимание формулу для периметра эллипса, как если бы эллипс был раздавлен. (Если малая ось в три раза меньше главной оси)

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

Периметр эллипса

Ответ 8

Существует действующий код MATLAB который можно найти здесь. Я повторяю это ниже, если эта связь когда-либо гаснет. Кредиты принадлежат оригинальному автору.

Этот код предполагает, что основная ось является отрезком от (x1, y1) до (x2, y2), а e является эксцентриситет эллипс.

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