[GIS] How to find the smallest-area ellipse circumscribing a convex polygon

algorithmgeometry

We are currently dealing with the problem of:
How to find an ellipse circumscribing a convex polygon?
Recent googling and lots of reading could not help me so much 🙁
enter image description here

Best Answer

An algorithm for your problem is described here. You may be interested by its implementation in CGAL.

PostGIS

If you are using the PostGIS stack, pgRouting's optional Driving Distance extension exposes an alpha shape implementation. I'm not sure if you can vary the alpha parameter, however.

Usage is below:

SELECT the_geom AS alpha_shape 
FROM 
  points_as_polygon(
    'SELECT id, ST_X(your_geom) AS x, ST_Y(your_geom) AS y FROM your_table');

Python

There are probably many Python modules you could use. I have heard good things about CGAL, a C++ computational geometry library. Python wrappers exist for parts of CGAL, including exposing CGAL's alpha shape implementation to Python.

Be aware that parts of CGAL are licensed under the QPL, which means that if you distribute your program, linked to CGAL, you may need to release it under the QPL. It is fine to keep your code proprietary if you do not redistribute your program code or binaries.

Geometry Center – Finding Center of Geometry of an Object

Every polygon has, at a minimum, four distinct "centers":

The barycenter of its vertices.
The barycenter of its edges.
Its barycenter as a polygon.
A GIS-specific "center" useful for labeling (usually calculated with undocumented proprietary methods).

(They may accidentally coincide in special cases, but for "generic" polygons they are distinct points.)

A "barycenter" in general is a "center of mass." The three types differ on where the mass is presumed located: it either is entirely on the vertices, spread uniformly on the edges, or spread uniformly throughout the polygon itself.

Simple methods exist to compute all three barycenters. One approach relies on the basic fact that the barycenter of the disjoint union of two masses is the total-mass-weighted average of the barycenters. From this we easily obtain the following:

The barycenter of two (equally weighted) vertices is their average. This is obtained by averaging their coordinates separately. Geometrically, it is the midpoint of the line segment joining the two vertices.
Inductively, the barycenter of n (equally weighted) vertices is obtained by averaging their coordinates separately.
The barycenter of a line segment is its midpoint. (This is clear by symmetry.)
The barycenter of a polyline is obtained by finding the midpoints of each line segment and then forming their weighted average using the segment lengths as weights.

For example, consider the "L" shape delineated by the points (0,0), (6,0), (6,12). There are two segments: one of length 6 with midpoint at ( (0+0)/2, (0+6)/2 ) = (3,0) and another of length 12 with midpoint at ( (6+6)/2, (0+12)/2 ) = (6,6). Their length-weighted average coordinates are therefore (x,y) with
```
x = (6*3 + 12*6) / (6+12) = 5,  y = (6*0 + 12*6) / (6+12) = 4.
```
This differs from the barycenter of the three vertices, which is ( (0+6+6)/3, (0+0+12)/3 ) = (4,4).

(Edit As another example, consider the figure in the question, which although square in shape, is represented as a pentagon determined by the sequence of points (0,0), (1/2,0), (1,0), (1,1), (0,1). The five sides have lengths 1/2, 1/2, 1, 1, 1 and midpoints (1/4,0), (3/4,0), (1,1/2), (1/2,1), and (0,1/2), respectively. Their weighted average therefore equals
```
[(1/2)*(1/4, 0) + (1/2)*(3/4, 0) + (1)*(1, 1/2) + (1)*(1/2, 1) + (1)*(0, 1/2)] / (1/2+1/2+1+1+1)
= (2/4, 2/4) = (0.5, 0.5)
```
as one would hope, even though the barycenter of the vertices alone (computed as in #2 above) is (0.5, 0.4).)
The barycenter of a polygon can be obtained by triangulation to decompose it into triangles. The barycenter of a triangle-qua-polygon coincides with the barycenter of its vertices. The area-weighted average of these barycenters is the polygon's barycenter. Triangle areas are readily computed in terms of their vertex coordinates (e.g., in terms of the wedge product of two of the sides). For an illustration of such area calculations, including how to exploit signed (positive or negative) areas, see the section on "Area" at my (old) course notes page.

(Edit Consider the polygon depicted in the question for example. We could triangulate it with triangles ((0,0), (1/2,0), (0,1)) on the left, ((0,1), (1/2,0), (1,1)) in the middle, and ((1,1), (1,0), (1/2,0)) on the right. Their areas are 1/4, 1/2, 1/4 respectively and their barycenters--obtained by averaging their vertices--are (1/6,1/3), (1/2,2/3), and (5/6,1/3), respectively. The area-weighted average of these barycenters equals
```
[(1/4)*(1/6,1/3) + (1/2)*(1/2,2/3) + (1/4)*(5/6,1/3)] / (1/4 + 1/2 + 1/4)
= (12/24, 6/12)
= (0.5, 0.5)
```
as it should, despite the presence of that fifth vertex along the bottom edge.)

It is evident that each of these methods is efficient: it requires just a single pass over the "spaghetti" representation of the polygon, using (fairly little) constant time at each step. Note that in all cases except the first (of pure vertices), more information than just a list of vertex coordinates is needed: you need to know the topology of the figure as well. In the "L" example, we needed to know that (0,0) was connected to (6,0) and not to (6,12), for instance.

These are all Euclidean concepts. They can be extended to the sphere (or ellipsoid) in several ways. A straightforward one views the features as a simplicial complex in three (Euclidean) dimensions, computes the appropriate barycenter, and then projects it outward from the center of the ellipsoid back to the surface. This requires no new concepts or formulas; you only have to work with a third (z) coordinate in addition to the first two coordinates. (Areas are still found using lengths of wedge products.)

Another generalization recognizes that the Euclidean metric--the square root of a sum of squares, according to Pythagoras--can be changed to other Lp metrics for p >= 1: you take the pth root of the sum of pth powers. Finding appropriate "barycenters" is no longer so simple, because the beautiful additive properties exploited above (barycenters are weighted averages of barycenters of simpler parts of a figure) no longer hold in general. Often, iterative approximate numerical solutions have to be obtained. They might not even be unique.

Additional centers can be defined for various purposes. Triangles have many different centers that can generalize (somewhat) to polygons: the center of the circumcircle, the center of (some) maximal incircle, the center of a minimum-area bounding ellipse, and others. Any set can be enclosed in various "hulls," such as the convex hull, and the centers of those hulls obtained.

Note that many of these "centers" are not necessarily located within the interior of a polygon. (Any reasonable center of a convex polygon will lie within its interior, though.)

This variety of approaches and solutions indicates one should be wary of a generic term like "center of geometry" or merely "center": it could be just about anything.

Best Answer

Related Solutions

Concave Hull Algorithm – Definition, Algorithms, and Practical Solutions for Concave Hull

PostGIS

Python

Geometry Center – Finding Center of Geometry of an Object

Related Question