用多边形逼近椭圆

Question

我正在处理地理信息，最近我需要绘制一个椭圆。为了与 OGC 约定兼容，我不能按原样使用椭圆；取而代之的是，我使用多边形来近似椭圆，方法是采用椭圆包含的多边形并使用任意多个点。

我用来为给定数量的点 N 生成椭圆的过程如下（使用 C# 和虚构的 Polygon 类）：

Polygon CreateEllipsePolygon(Coordinate center, double radiusX, double radiusY, int numberOfPoints)
{
    Polygon result = new Polygon();
    for (int i=0;i<numberOfPoints;i++)
    {
        double percentDone = ((double)i)/((double)numberOfPoints);
        double currentEllipseAngle = percentDone * 2 * Math.PI;
        Point newPoint = CalculatePointOnEllipseForAngle(currentEllipseAngle, center, radiusX, radiusY);
        result.Add(newPoint);
    }
    return result;
}

到目前为止，这已经为我服务了很长时间，但我注意到它有一个问题：如果我的椭圆是“粗壮”的，也就是说，radiusX 远大于 radiusY，则顶部的点数椭圆与椭圆左边的点数相同。

这是对积分的浪费！在椭圆的上部添加一个点几乎不会影响我的多边形逼近的精度，但在椭圆的左侧添加一个点会产生重大影响。

我真正想要的是一种用多边形近似椭圆的更好算法。我需要从这个算法中得到什么：

• 它必须接受点数作为参数；可以接受每个象限中的点数（我可以迭代地在“有问题”的地方添加点，但我需要很好地控制我正在使用的点数）
• 必须以椭圆为界
• 它必须包含椭圆中心的正上方、正下方、向左和向右的点
• 它的面积应该尽可能接近椭圆的面积，当然优先考虑给定点数的最优值（见 Jaan 的回答 - 显然这个解决方案已经是最优的了）
• 多边形的最小内角最大

我想到的是找到一个多边形，其中每两条线之间的角度始终相同 - 但我不仅不知道如何生成这样的多边形，我什至不确定存在，即使我取消了限制！

有人知道如何找到这样的多边形吗？

Answer 1

我假设在 OP 的问题中，CalculatePointOnEllipseForAngle 返回一个坐标如下的点。

newPoint.x = radiusX*cos(currentEllipseAngle) + center.x
newPoint.y = radiusY*sin(currentEllipseAngle) + center.y

那么，如果目标是最小化椭圆和内接多边形的面积差异（即找到一个面积最大的内接多边形），则OP的原始解决方案已经是最优的了。见Ivan Niven, "Maxima and Minima Without Calculus", Theorem 7.3b。 （有无限多个最优解：可以通过在上面的公式中向currentEllipseAngle 添加任意常数来获得另一个具有相同面积的多边形；这些是唯一的最优解。证明思想很简单：第一个证明这些是圆形情况下的最佳解决方案，即如果radiusX=radiusY；其次观察到，在将圆形转换为椭圆的线性变换下，例如将x坐标乘以某个常数的变换，所有面积乘以一个常数，因此圆的最大面积内接多边形转换为椭圆的最大面积内接多边形。）

人们也可以考虑其他目标，如其他帖子中所建议的那样：例如最大化多边形的最小角度或最小化多边形和椭圆边界之间的Hausdorff distance。 （例如，Ramer-Douglas-Peucker algorithm 是一种近似解决后一个问题的启发式方法。不像在通常的 Ramer-Douglas-Peucker 实现中那样近似多边形曲线，我们近似一个椭圆，但可以设计一个公式来查找一个椭圆是离线段最远的点。）关于这些目标，OP的解决方案通常不是最优的，我不知道找到一个精确的解决方案公式是否可行。但是 OP 的解决方案并不像 OP 的图片显示的那么糟糕：似乎 OP 的图片不是使用该算法生成的，因为它在椭圆的更尖锐弯曲部分中的点比该算法生成的要少。

Answer 2

finding a polygon in which the angle between every two lines is
always the same

是的，这是可能的。我们想找到（第一个）椭圆象限的这样的点，这些点的切线角度形成等距（相同的角度差）序列。不难找到切点

x=a*Cos(fi)
y=b*Sin(Fi)

derivatives
dx=-a*Sin(Fi), dy=b*Cos(Fi)
y'=dy/dx=-b/a*Cos(Fi)/Sin(Fi)=-b/a*Ctg(Fi) 

导数y'描述切线，这个切线有角系数

k=b/a*Cotangent(Fi)=Tg(Theta)
Fi = ArcCotangent(a/b*Tg(Theta)) = Pi/2-ArcTan(a/b*Tg(Theta)) 

由于relation for complementary angles

其中 Fi 从 0 变化到 Pi/2，而 Theta - 从 Pi/2 变化到 0。
因此，每个象限查找 N + 1 个点（包括极值点）的代码可能看起来像（这是生成附加图片的 Delphi 代码）

  for i := 0 to N - 1 do begin
    Theta := Pi/2 * i /  N;
    Fi :=  Pi/2 - ArcTan(Tan(Theta) * a/b);
    x := CenterX + Round(a * Cos(Fi));
    y := CenterY + Round(b * Sin(Fi));
  end;
  // I've removed Nth point calculation, that involves indefinite Tan(Pi/2) 
  // It would better to assign known value 0 to Fi in this point

完美角多边形的草图：

Answer 3

这是我用过的迭代算法。

我没有寻找理论上最优的解决方案，但它对我来说效果很好。

请注意，此算法将多边形素数对椭圆的最大误差作为输入，而不是您希望的点数。

public static class EllipsePolygonCreator
{
#region Public static methods

public static IEnumerable<Coordinate> CreateEllipsePoints(
  double maxAngleErrorRadians,
  double width,
  double height)
{
  IEnumerable<double> thetas = CreateEllipseThetas(maxAngleErrorRadians, width, height);
  return thetas.Select(theta => GetPointOnEllipse(theta, width, height));
}

#endregion

#region Private methods

private static IEnumerable<double> CreateEllipseThetas(
  double maxAngleErrorRadians,
  double width,
  double height)
{
  double firstQuarterStart = 0;
  double firstQuarterEnd = Math.PI / 2;
  double startPrimeAngle = Math.PI / 2;
  double endPrimeAngle = 0;

  double[] thetasFirstQuarter = RecursiveCreateEllipsePoints(
    firstQuarterStart,
    firstQuarterEnd,
    maxAngleErrorRadians,
    width / height,
    startPrimeAngle,
    endPrimeAngle).ToArray();

  double[] thetasSecondQuarter = new double[thetasFirstQuarter.Length];
  for (int i = 0; i < thetasFirstQuarter.Length; ++i)
  {
    thetasSecondQuarter[i] = Math.PI - thetasFirstQuarter[thetasFirstQuarter.Length - i - 1];
  }

  IEnumerable<double> thetasFirstHalf = thetasFirstQuarter.Concat(thetasSecondQuarter);
  IEnumerable<double> thetasSecondHalf = thetasFirstHalf.Select(theta => theta + Math.PI);
  IEnumerable<double> thetas = thetasFirstHalf.Concat(thetasSecondHalf);
  return thetas;
}

private static IEnumerable<double> RecursiveCreateEllipsePoints(
  double startTheta,
  double endTheta,
  double maxAngleError,
  double widthHeightRatio,
  double startPrimeAngle,
  double endPrimeAngle)
{
  double yDelta = Math.Sin(endTheta) - Math.Sin(startTheta);
  double xDelta = Math.Cos(startTheta) - Math.Cos(endTheta);
  double averageAngle = Math.Atan2(yDelta, xDelta * widthHeightRatio);

  if (Math.Abs(averageAngle - startPrimeAngle) < maxAngleError &&
      Math.Abs(averageAngle - endPrimeAngle) < maxAngleError)
  {
    return new double[] { endTheta };
  }

  double middleTheta = (startTheta + endTheta) / 2;
  double middlePrimeAngle = GetPrimeAngle(middleTheta, widthHeightRatio);
  IEnumerable<double> firstPoints = RecursiveCreateEllipsePoints(
    startTheta,
    middleTheta,
    maxAngleError,
    widthHeightRatio,
    startPrimeAngle,
    middlePrimeAngle);
  IEnumerable<double> lastPoints = RecursiveCreateEllipsePoints(
    middleTheta,
    endTheta,
    maxAngleError,
    widthHeightRatio,
    middlePrimeAngle,
    endPrimeAngle);

  return firstPoints.Concat(lastPoints);
}

private static double GetPrimeAngle(double theta, double widthHeightRatio)
{
  return Math.Atan(1 / (Math.Tan(theta) * widthHeightRatio)); // Prime of an ellipse
}

private static Coordinate GetPointOnEllipse(double theta, double width, double height)
{
  double x = width * Math.Cos(theta);
  double y = height * Math.Sin(theta);
  return new Coordinate(x, y);
}

#endregion
}

Answer 4

为闭合轮廓（如椭圆）实现自适应离散化的一种方法是反向运行Ramer–Douglas–Peucker algorithm：

1. Start with a coarse description of the contour C, in this case 4 
   points located at the left, right, top and bottom of the ellipse.
2. Push the initial 4 edges onto a queue Q.

while (N < Nmax && Q not empty)

3. Pop an edge [pi,pj] <- Q, where pi,pj are the endpoints.
4. Project a midpoint pk onto the contour C. (I expect that 
   simply bisecting the theta endpoint values will suffice
   for an ellipse).
5. Calculate distance D between point pk and edge [pi,pj].

    if (D > TOL)

6.      Replace edge [pi,pj] with sub-edges [pi,pk], [pk,pj].
7.      Push new edges onto Q.
8.      N = N+1

    endif

endwhile

该算法迭代地细化轮廓C 的初始离散化，将高曲率区域中的点聚类。它在满足(i) 用户定义的容错TOL 或(ii) 使用最大允许点数Nmax 时终止。

我确信可以找到一种专门针对椭圆的情况进行优化的替代方法，但我认为这种方法的通用性非常方便。

Answer 5

我建议你切换到极坐标：

极坐标椭圆是：

x(t) = XRadius * cos(t)
y(t) = YRadius * sin(t)

为0 <= t <= 2*pi

Xradius >> YRadius（或Yradius >> Yradius）出现问题

除了使用 numberOfPoints 之外，您还可以使用角度数组，显然它们并不完全相同。 IE。用 36 分平均划分每个扇区，您将得到 angle = 2*pi*n / 36 radiants。 当您在这两个值的“邻域”中获得 n = 0（或 36）或 n = 18 时，近似方法不能很好地工作，因为椭圆扇区与用于近似它的三角形有很大不同。您可以减小此点周围的扇区大小，从而提高精度。而不仅仅是增加点的数量，这也会增加其他不需要的区域中的段。角度的顺序应该变成（以度为单位）：

angles_array = [5,10,10,10,10.....,5,5,....10,10,...5]

前 5 度。序列是 t = 0，第二个是 t = pi，最后一个是大约 2*pi。