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:
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.
Remove all duplicate vertices, including the end vertex if it is the same as the start vertex.
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.
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.
After each run is identified, append its end vertex to the new polygon.
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.
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.