如何生成沿椭圆周长均匀分布的一组点？
Eri*_*ric 13 language-agnostic math geometry
如果我想生成一堆围绕圆形均匀分布的点,我可以这样做(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))
]
是否有一种简单的方法可以在椭圆周围生成等间距的点？
小智 12
这是一个古老的线程,但因为我寻求沿或椭圆,创造均匀间隔点的相同任务不是能够找到一个实现,我提供了实现霍华德的伪代码这个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;
}
}
（更新：以反映新包装）。
可以在Python 包的 numeric 分支FlyingCircus-Numeric中找到 Python 的此问题的有效解决方案FlyingCircus。
免责声明：我是它们的主要作者。
简而言之，（简化的）代码看起来（哪里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()
可能的（数值）计算如下：
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
这是一个简单的数值积分方案。如果您需要更高的准确性，您也可以使用任何其他集成方法。
我确定这个线程现在已经死了很久了，但我刚刚遇到了这个问题，这是最接近解决方案的。
我从戴夫在这里的回答开始，但我注意到它并没有真正回答海报的问题。它不是按弧长等分椭圆，而是按角度。
无论如何，我对他的（很棒的）工作做了一些调整，让椭圆被弧长等分（这次是用 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;
}
| 归档时间： |
|
| 查看次数： |
8311 次 |
| 最近记录： |