# Original polygon is completely contained by new after simplifying

Is there a polygon simplification method that can create a polygon, such that the original is completely contained by the new one?

The aim is to create the new polygon from as few points as possible while maintaining a specified maximum tolerance. Until now I have been using a method that creates a buffer to a distance of specified tolerance and then uses a regular simplification method but that seems to be very inefficient as it sometimes creates an unnecessary offset. The original does not intersect itself.

My desired output looks as following

Simplify your original polygon first. Than use the offset curves tool (from advanced digitizing toolbox - activate it if not visible) to enlarge the simplified polygon. Activate the snapping-mode and make sure that the enlarged polygon snaps to the point on your original polygon that is farthest away from your simplified one (to avoid unnecessary offset).

If offset curves produces round angles configure the parameters in QGIS options dialogue

See screenshot:

I've had to solve this problem recently in my commercial cartography library. I won't name it because we're not suppose to advertise here.

The reason for needing the solution is to have a quick way to see if a point is in a polygon, at the risk of some false positives.

The function takes a polygon consisting of an ordered sequence of 2D vertices. Integer arithmetic can be used as long as care is taken to avoid overflow when calculating cross-products. There is one parameter: MaxNewEdgeLength, which is the maximum length of any new edge that is inserted. If you set that number to something very large, a standard convex hull is produced.

The algorithm is:

1. Find out whether the polygon is clockwise or anti-clockwise using the standard method: find the bottom point (taking the leftmost of any duplicates) and see whether the partial area, determined by the cross product of the vectors before and after it, is positive or negative. Positive means anti-clockwise. Proviso: you may want to flip the result for inner polygons if you're dealing with a poly-polygon.

2. Remove all duplicate vertices, including the end vertex if it is the same as the start vertex.

3. Label every vertex as concave, flat or convex using the cross-product method described above, in combination with whether the polygon is clockwise or anticlockwise. For an anticlockwise polygon, a negative cross-product means convex and a positive one means concave; it's the reverse for clockwise polygons. A cross product of zero means a flat vertex; you can delete all flat vertices if you like, to simplify what follows.

4. Optionally (it doesn't make much difference) skip forward to the first convex vertex. Then traverse the polygon, finding runs of these types:

• A flat run is a run of flat vertices, plus the non-flat vertex before and after it.
• A concave run is a run of concave vertices, plus the non-concave vertex before and after it, BUT a concave run is terminated if extending it would make the distance between start and end vertices greater than MaxNewEdgeLength.
• A convex run is a run consisting of a convex vertex and the vertex before it.
1. After each run is identified, append its end vertex to the new polygon.

2. Has the new polygon fewer vertices than the original one? If yes, try another iteration by going back to step 3 and using the new polygon as input. You may want to limit the number of iterations: I use an arbitrary limit of 32, which works well with OpenStreetMap cartographic data.

3. Return the new polygon.

Now some caveats: (i) the algorithm is not optimal, in that it is clear that there are better solutions. That's because it takes the largest bite possible, limited by MaxNewEdgeLength, when finding a concave run, and quite often it would be better to split a large concave run in a different place. (ii) I'm morally certain that it could be coded up in a faster way. (iii) It will do nothing for convex polygons. These almost never occur in the data I'm working with, which is administrative boundaries.

Here are some images showing the OpenStreetMap outline of the city of Bath, England, reduced from 1454 points to 96, using a MaxNewEdgeLength of 1km. The original points are shown as red (convex) and green (concave). The simplified outline is drawn in blue.

Code to do exactly this is now available in the JTS Topology Suite. It is in the `PolygonHull` class. It supports both Outer and Inner hulls. The control parameterization is either fraction of vertices in result or ratio of change in area.

See the PR.