17

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.

Illustration

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?

Gilthans
  • 1,656
  • 1
  • 18
  • 23
  • 1
    I suppose one way to approach this problem is to consider the gradient (derivative) of the boundary of you polygon. Consider a point on your polygon, and its two direct neighbors. You can compute the angle that is formed by the two edges emanating from the point to it's neighbors. The more points you have, that angle will tend to 180 degrees. When you have few points, the points that are "unnecessary" will have large angles (close to 180), the "needs more" will have small angles. So add more points there, remove the unneeded once. Hope that helps. – Nicko Po Mar 27 '14 at 17:44
  • 1
    What you suggest is similar to what MaMazav offered below, which like I said is a good fallback, but I'd prefer something with more guarantees – Gilthans Mar 29 '14 at 13:20
  • The picture has been drawn manually, not by the algorithm, has it? The picture exaggregates the situation worse than it is, if `CalculatePointOnEllipseForAngle` is implemented as `new Point (radiusX*cos(currentEllipseAngle) + center.x, radiusY*sin(currentEllipseAngle) + center.y)`. – Jaan Mar 26 '15 at 21:29
  • I know this is quite old at this point, but I spent hours searching for a solution and the MBo response was close for me, but didn't quite work when I implemented the code example as shown. I found this example which I translated into JS and worked great! https://www.geeksforgeeks.org/how-to-discretize-an-ellipse-or-circle-to-a-polygon-using-c-graphics/ specifically the code discretizeEllipse, it works great!, as you can see it is VERY similar to MBo's response. – Krum110487 Jun 01 '22 at 17:53

5 Answers5

14
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

enter image description here

Sketch for perfect-angle polygon:

enter image description here

MBo
  • 77,366
  • 5
  • 53
  • 86
  • This looks awesome, +1 for the code sample and the attached image. It evens fits perfectly to my current design. I'll try this and accept the answer if it works :) You mentiond Fi and how it is calculated, but I didn't understand why it has this value, it seems an almost random manipulation of the angle. Can you elaborate on how you chose Fi? – Gilthans Mar 28 '14 at 16:03
  • The part where you transformed ArcCot wasn't the problem, the problem is more about how you found k, and what does it mean geometrically. It seems like you've divided y by x (and I don't know why), but then you'd have gotten Tg, not Cotg. Also, why would that be equal to Tg(theta)? – Gilthans Mar 28 '14 at 22:21
  • I've tried implementing the algorithm. It gave a pretty result in that the points were denser around the edges, but after some calculation and observation I noticed the polygon's internal angle were not equal, with a variance of up to 3 degrees, which is pretty significant. It seems like this algorithm is TOO biased towards the edges - I've noticed now that the ellipse's vertical sides have unnecessary points, while the ellipse's horizontal sides are missing a few. I'd try and tweak the algorithm, if I understood it :/ – Gilthans Mar 28 '14 at 22:23
  • Polygon angles will be exactly equal, if we build external tangents to ellipse in calculated points, and find their pairwise intersections. Picture added. – MBo Mar 29 '14 at 02:12
  • I can't attach a picture or code, but here's the result I'm getting using this algorithm: http://i.imgur.com/tZxss3c.png I think the problem has to do with the fact that it takes into account every quadrant in itself, and so the bottom point for example gets twice as much acuteness than it should. It's cool in the sense that it is a step in the direction I want, but the result is far from optimal and doesn't look like it fits my bill – Gilthans Mar 29 '14 at 13:18
  • Of course, extremal points have to be used once, not twice. Your picture - it is in the nature of things for inscribed polygon with small amount of edges. You can add more points or try escribed polygon as bottom picture shows – MBo Mar 29 '14 at 13:58
  • I don't see how this is a property of inscribed polygons... There is a clear optimization without adding any new points, by 'moving' the unnecessary point to where it is needed, suggesting that this solution is not optimal. Also, it doesn't seem to have equal angles, as confirmed both by eye sight and calculation. In your example there are many points so it looks pretty good, but the algorithm is intended to give good result for as little points as possible – Gilthans Mar 29 '14 at 15:04
  • Equal angles - for tangents (edges od escribed polygon), not for edges of inscribed polygon. I believed that inscribed polygon is quite good approximation with a reasonable amount of points. For perfectly equal angles - use escribed polygon with edge centers in given points, or inscribed polygon calculated with another method (constant edge-edge angle) – MBo Mar 29 '14 at 15:24
  • The problem with this solution (which is very good) is that it gives an optimal solution for an escribed polygon, whereas I need the inscribed polygon. Otherwise this is a great answer. – Gilthans May 28 '14 at 21:03
  • 1
    I got much more pleasing results by changing the Fi := Pi/2 - ArcTan(Tan(Theta) * a/b); line to Fi := Pi/2 - ArcTan(Tan(Theta) * Sqrt(a/b)); – mheyman Aug 20 '15 at 14:30
  • and how big should N be (I understand that question says "given number of points", but I'm interested in more general case) – Superior May 27 '22 at 14:02
  • @Superior Depends on ellipse size and desired quality. I think that about 10 points per quadrant give visible kinks, 100 points should be good looking – MBo May 27 '22 at 14:09
3

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.

Darren Engwirda
  • 6,915
  • 4
  • 26
  • 42
  • This is interesting but seems far from ideal. Unless I select TOL very well, the algorithm will randomly select edges and find their midpoint. My intuition says this will return the same result I've had so far unless I somehow find a good TOL value (which I don't know) – Gilthans Mar 28 '14 at 22:18
  • The result would be a little better and TOL would have a simple meaning (the distance between the farthest point of the ellipse from the polygon, e.g. 1 pixel or 5 meters or whatever your requirements for the allowed error are), if instead of projecting "a midpoint pk onto the contour C" in point 4, one would find the point on the arc of ellipse that is farthest from the polygon's edge (with some calculus, one can derive a formula for that). Then it could be called the Douglas–Peucker algorithm itself, not "Douglas–Peucker algorithm in reverse". – Jaan Aug 19 '15 at 11:08
3

I assume that in the OP's question, CalculatePointOnEllipseForAngle returns a point whose coordinates are as follows.

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

Then, if the goal is to minimize the difference of the areas of the ellipse and the inscribed polygon (i.e., to find an inscribed polygon with maximal area), the OP's original solution is already an optimal one. See Ivan Niven, "Maxima and Minima Without Calculus", Theorem 7.3b. (There are infinitely many optimal solutions: one can get another polygon with the same area by adding an arbitrary constant to currentEllipseAngle in the formulae above; these are the only optimal solutions. The proof idea is quite simple: first one proves that these are the optimal solutions in case of a circle, i.e. if radiusX=radiusY; secondly one observes that under a linear transformation that transforms a circle into our ellipse, e.g. a transformation of multiplying the x-coordinate by some constant, all areas are multiplied by a constant and therefore a maximal-area inscribed polygon of the circle is transformed into a maximal-area inscribed polygon of the ellipse.)

One may also regard other goals, as suggested in the other posts: e.g. maximizing the minimal angle of the polygon or minimizing the Hausdorff distance between the boundaries of the polygon and ellipse. (E.g. the Ramer-Douglas-Peucker algorithm is a heuristic to approximately solve the latter problem. Instead of approximating a polygonal curve, as in the usual Ramer-Douglas-Peucker implementation, we approximate an ellipse, but it is possible to devise a formula for finding on an ellipse arc the farthest point from a line segment.) With respect to these goals, the OP's solution would usually not be optimal and I don't know if finding an exact solution formula is feasible at all. But the OP's solution is not as bad as the OP's picture shows: it seems that the OP's picture has not been produced using this algorithm, as it has less points in the more sharply curved parts of the ellipse than this algorithm produces.

Jaan
  • 2,220
  • 20
  • 17
  • +1 for referencing a proof (albeit one I don't have access to). The goal indeed has to change in this case - I think maximizing the minimal angle is a good goal. While the picture indeed wasn't generated (I don't have access to the original code for various reasons), the problem grows increasingly worse the more 'stocky' the ellipse is (IE larger difference between x-radius and y-radius) – Gilthans Aug 19 '15 at 18:59
  • @Gilthans I added an explanation of the proof idea for those who don't have access to the book link. – Jaan Aug 23 '15 at 14:31
1

I suggest you switch to polar coordinates:

Ellipse in polar coord is:

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

for 0 <= t <= 2*pi

The problems arise when Xradius >> YRadius (or Yradius >> Yradius)

Instead of using numberOfPoints you can use an array of angles obviously not all identical. I.e. with 36 points and dividing equally you get angle = 2*pi*n / 36 radiants for each sector. When you get around n = 0 (or 36) or n = 18 in a "neighborhood" of these 2 values the approx method doesn't works well cause the ellipse sector is significantly different from the triangle used to approximate it. You can decrease the sector size around this points thus increasing precision. Instead of just increasing the number of points that would also increase segments in other unneeded areas. The sequence of angles should become something like (in degrees ):

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

The first 5 deg. sequence is for t = 0 the second for t = pi, and again the last is around 2*pi.

EWit
  • 1,954
  • 13
  • 22
  • 19
Max
  • 11
  • 1
1

Here is an iterative algorithm I've used.

I didn't look for theoretically-optimal solution, but it works quit well for me.

Notice that this algorithm gets as an input the maximal error of the prime of the polygon agains the ellipse, and not the number of points as you wish.

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
}
MaMazav
  • 1,773
  • 3
  • 19
  • 33