Approximating an ellipse with a polygon

Question

I am working with geographic information, and recently I needed to draw an ellipse. For compatibility with the OGC convention, I cannot use the ellipse as it is; instead, I use an approximation of the ellipse using a polygon, by taking a polygon which is contained by the ellipse and using arbitrarily many points.

The process I used to generate the ellipse for a given number of point N is the following (using C# and a fictional Polygon class):

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

This has served me quite while so far, but I've noticed a problem with it: if my ellipse is 'stocky', that is, radiusX is much larger than radiusY, the number of points on the top part of the ellipse is the same as the number of points on the left part of the ellipse.

That is a wasteful use of points! Adding a point on the upper part of the ellipse would hardly affect the precision of my polygon approximation, but adding a point to the left part of the ellipse can have a major effect.

What I'd really like, is a better algorithm to approximate the ellipse with a polygon. What I need from this algorithm:

• It must accept the number of points as a parameter; it's OK to accept the number of points in every quadrant (I could iteratively add points in the 'problematic' places, but I need good control on how many points I'm using)
• It must be bounded by the ellipse
• It must contain the points straight above, straight below, straight to the left and straight to the right of the ellipse's center
• Its area should be as close as possible to the area of the ellipse, with preference to optimal for the given number of points of course (See Jaan's answer - appearantly this solution is already optimal)
• The minimal internal angle in the polygon is maximal

What I've had in mind is finding a polygon in which the angle between every two lines is always the same - but not only I couldn't find out how to produce such a polygon, I'm not even sure one exists, even if I remove the restrictions!

Does anybody have an idea about how I can find such a polygon?

Answer 1

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

Yes, it is possible. We want to find such points of (the first) ellipse quadrant, that angles of tangents in these points form equidistant (the same angle difference) sequence. It is not hard to find that tangent in point

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) 

Derivative y' describes tangent, this tangent has angular coefficient

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

due to relation for complementary angles

where Fi varies from 0 to Pi/2, and Theta - from Pi/2 to 0.
So code for finding N + 1 points (including extremal ones) per quadrant may look like (this is Delphi code producing attached picture)

  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

Sketch for perfect-angle polygon:

Answer 2

One way to achieve adaptive discretisations for closed contours (like ellipses) is to run the Ramer–Douglas–Peucker algorithm in reverse:

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

This algorithm iteratively refines an initial discretisation of the contour C, clustering points in areas of high curvature. It terminates when, either (i) a user defined error tolerance TOL is satisfied, or (ii) the maximum allowable number of points Nmax is used.

I'm sure that it's possible to find an alternative that's optimised specifically for the case of an ellipse, but the generality of this method is, I think, pretty handy.