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?

4
  • 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, 2014 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, 2014 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, 2015 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! geeksforgeeks.org/… specifically the code discretizeEllipse, it works great!, as you can see it is VERY similar to MBo's response.
    – Krum110487
    Jun 1, 2022 at 17:53

5 Answers 5

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

12
  • 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, 2014 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, 2014 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, 2014 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, 2014 at 2:12
  • 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, 2015 at 14:30
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.

2
  • 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, 2014 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, 2015 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.

2
  • +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, 2015 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, 2015 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.

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
}
1
  • This kind of solution will be my fallback, but I'd prefer something with more guarantees
    – Gilthans
    Mar 29, 2014 at 9:10

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.