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

Если я хочу создать группу точек, равномерно распределенных по кругу, я могу сделать это (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))
]

Существует ли простой способ создания точек на одинаковом расстоянии вокруг эллипса?

Ответы [ 8 ]

12 голосов
/ 11 декабря 2013

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

}
10 голосов
/ 07 августа 2011

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

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

4 голосов
/ 07 августа 2011

Возможный (числовой) расчет может выглядеть следующим образом:

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

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

2 голосов
/ 28 августа 2018

Эффективное решение этой проблемы для 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()

angles_in_ellipse

2 голосов
/ 15 июня 2015

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

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

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

    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;
    }
0 голосов
/ 06 июля 2019

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

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

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

0 голосов
/ 02 ноября 2016

Работающий код 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
0 голосов
/ 10 августа 2016

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

Я добавляю его в 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()
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...