67

Having a set of (2D) points from a GIS file (a city map), I need to generate the polygon that defines the 'contour' for that map (its boundary). Its input parameters would be the points set and a 'maximum edge length'. It would then output the corresponding (probably non-convex) polygon.

The best solution I found so far was to generate the Delaunay triangles and then remove the external edges that are longer than the maximum edge length. After all the external edges are shorter than that, I simply remove the internal edges and get the polygon I want. The problem is, this is very time-consuming and I'm wondering if there's a better way.

StayOnTarget
  • 11,743
  • 10
  • 52
  • 81
Fabio Ceconello
  • 15,819
  • 5
  • 38
  • 51
  • You say you have a gis file - are you not using a GIS mapping application / software? I know that ArcGIS server will happily consume any number of points and draw up a map overlain with the resulting polygon. – Ian Sep 17 '08 at 14:21
  • Yes, I have a GIS file but I need to write the algorithm (in C or C++, probably), this is to be placed inside an already existing program and it shouldn't use external tools (like ArcGIS) to do that, it needs to be self-contained. – Fabio Ceconello Sep 17 '08 at 14:26
  • 1
    Actually, I don't think ArcGIS has a built-in algorithm to do what he wants. ArcGIS has the ability to do convex hulls, but concave ones are considerably more complicated. – Chris Upchurch Sep 17 '08 at 14:50
  • 2
    Could you define your problem more precisely? :) With 5 points: 4 corners of a square, and its centre. What would your boundary be? If your maximum edge length allowed the centre, it's completely arbitrary as to which of the 4 edges of the square you would 'bend in' to include the middle point. – Luke Halliwell Sep 17 '08 at 17:25
  • In the example you gave, the answer would always be the square. As a general rule, given each point, and having the distances for the 2 nearest other points, you can assume that the maximum edge length would never be smaller than the second one - so there should never be 'orphaned' points. – Fabio Ceconello Sep 17 '08 at 23:45
  • @Fabio Ceconello - Did you ever implement an algorithm for this? – stevehipwell Jan 29 '10 at 12:07
  • 1
    Yes, I did, the way I mentioned in the question (with the Delaunay triangles). Later I did some work on the alpha shapes concept pointed by nsanders below, but before I made it work the issue got its priority lowered and I moved to another one. – Fabio Ceconello Feb 07 '10 at 00:57
  • 2
    Delaunay has correct complexity (O(n log n)). I doubt you can do much better asymptotically. – Alexandre C. Aug 28 '11 at 10:27
  • You've accepted an answer after just 9 years! – peterh Mar 24 '17 at 16:16
  • @peterh I had accepted an answer right away, but it was removed later for some reason, so when I noticed that I accepted another - that would be my "second best". – Fabio Ceconello Mar 24 '17 at 16:31
  • @peterh I'm of the opinion no answer should be accepted until someone provides proper code. If Fabio posts even a snippet his code, I would upvote his over any of these answers and downvote all the rest to move it up the list – Nathan majicvr.com Feb 01 '19 at 16:28
  • Also, I'm confused about your solution. Did you have to manually remove edges, or was your solution automatic? (ie. code solved the whole thing end-to-end) At the very least, it sounds like you had to pick the "maximum edge length" by trial and error, rather than in any algorithmic way – Nathan majicvr.com Feb 01 '19 at 16:40
  • @peterh thinking about it I guess you're right, I'll leave it without an accepted answer. About the snippet, I'll have to revisit some old code and I don't even remember how small (or big) such snippet would be and how much I'll have to edit to make it understandable. If feasible, I'll post it. – Fabio Ceconello Feb 01 '19 at 18:26
  • 1
    @frank It was not manual, it was part of an automated tool to process maps for a GPS nav app. It was indeed arbitrary in my specific case - the points were street corners, and the resulting polygon would be the contour for a city. I used an arbitrary value that would give me a detailed enough polygon that wouldn't be too heavy to render. I think it has to be this way, you have to determine the maximum length according to the needs of your application - I don't see how you could calculate it automatically beforehand. – Fabio Ceconello Feb 01 '19 at 18:36
  • @FabioCeconello Noo! If the answer answers the question, then accept it! I was so happy that I found the longest accepted question, I didn't want you unaccept it! – peterh Feb 01 '19 at 23:15

10 Answers10

11

One of the former students in our lab used some applicable techniques for his PhD thesis. I believe one of them is called "alpha shapes" and is referenced in the following paper:

http://www.cis.rit.edu/people/faculty/kerekes/pdfs/AIPR_2007_Gurram.pdf

That paper gives some further references you can follow.

nsanders
  • 12,250
  • 2
  • 40
  • 47
4

This paper discusses the Efficient generation of simple polygons for characterizing the shape of a set of points in the plane and provides the algorithm. There's also a Java applet utilizing the same algorithm here.

Amirali
  • 116
  • 7
3

The guys here claim to have developed a k nearest neighbors approach to determining the concave hull of a set of points which behaves "almost linearly on the number of points". Sadly their paper seems to be very well guarded and you'll have to ask them for it.

Here's a good set of references that includes the above and might lead you to find a better approach.

Vinko Vrsalovic
  • 330,807
  • 53
  • 334
  • 373
  • 2
    It looks like this is the paper in question: http://repositorium.sdum.uminho.pt/bitstream/1822/6429/1/ConcaveHull_ACM_MYS.pdf The idea is exceptionally clever and simple--it's just a straightforward modification of the Graham scan technique for convex nulls, as far as I understand it. No need for finding a Delaunay triangulation, etc. – Zach Conn Mar 28 '14 at 21:33
  • 1
    Said paper name is "Concave hull: a k-nearest neighbours approach for the computation of the region occupied by a set of points" It's available here: http://repositorium.sdum.uminho.pt/handle/1822/6429 – deshu Dec 24 '19 at 01:14
2

The answer may still be interesting for somebody else: One may apply a variation of the marching square algorithm, applied (1) within the concave hull, and (2) then on (e.g. 3) different scales that my depend on the average density of points. The scales need to be int multiples of each other, such you build a grid you can use for efficient sampling. This allows to quickly find empty samples=squares, samples that are completely within a "cluster/cloud" of points, and those, which are in between. The latter category then can be used to determine easily the poly-line that represents a part of the concave hull.

Everything is linear in this approach, no triangulation is needed, it does not use alpha shapes and it is different from the commercial/patented offering as described here ( http://www.concavehull.com/ )

monnoo
  • 301
  • 3
  • 3
1

You can do it in QGIS with this plug in; https://github.com/detlevn/QGIS-ConcaveHull-Plugin

Depending on how you need it to interact with your data, probably worth checking out how it was done here.

Cameron
  • 169
  • 6
1

As a wildly adopted reference, PostGIS starts with a convexhull and then caves it in, you can see it here.

https://github.com/postgis/postgis/blob/380583da73227ca1a52da0e0b3413b92ae69af9d/postgis/postgis.sql.in#L5819

Evan Carroll
  • 78,363
  • 46
  • 261
  • 468
1

The Bing Maps V8 interactive SDK has a concave hull option within the advanced shape operations.

https://www.bing.com/mapspreview/sdkrelease/mapcontrol/isdk/advancedshapeoperations?toWww=1&redig=D53FACBB1A00423195C53D841EA0D14E#JS

Within ArcGIS 10.5.1, the 3D Analyst extension has a Minimum Bounding Volume tool with the geometry types of concave hull, sphere, envelope, or convex hull. It can be used at any license level.

There is a concave hull algorithm here: https://github.com/mapbox/concaveman

Jaybird64
  • 19
  • 3
1

A quick approximate solution (also useful for convex hulls) is to find the north and south bounds for each small element east-west.

Based on how much detail you want, create a fixed sized array of upper/lower bounds. For each point calculate which E-W column it is in and then update the upper/lower bounds for that column. After you processed all the points you can interpolate the upper/lower points for those columns that missed.

It's also worth doing a quick check beforehand for very long thin shapes and deciding wether to bin NS or Ew.

Martin Beckett
  • 94,801
  • 28
  • 188
  • 263
1

Good question! I haven't tried this out at all, but my first shot would be this iterative method:

  1. Create a set N ("not contained"), and add all points in your set to N.
  2. Pick 3 points from N at random to form an initial polygon P. Remove them from N.
  3. Use some point-in-polygon algorithm and look at points in N. For each point in N, if it is now contained by P, remove it from N. As soon as you find a point in N that is still not contained in P, continue to step 4. If N becomes empty, you're done.
  4. Call the point you found A. Find the line in P closest to A, and add A in the middle of it.
  5. Go back to step 3

I think it would work as long as it performs well enough — a good heuristic for your initial 3 points might help.

Good luck!

mattdm
  • 2,082
  • 25
  • 39
Rob Dickerson
  • 1,459
  • 9
  • 6
1

A simple solution is to walk around the edge of the polygon. Given a current edge om the boundary connecting points P0 and P1, the next point on the boundary P2 will be the point with the smallest possible A, where

H01 = bearing from P0 to P1
H12 = bearing from P1 to P2
A = fmod( H12-H01+360, 360 )
|P2-P1| <= MaxEdgeLength

Then you set

P0 <- P1
P1 <- P2

and repeat until you get back where you started.

This is still O(N^2) so you'll want to sort your pointlist a little. You can limit the set of points you need to consider at each iteration if you sort points on, say, their bearing from the city's centroid.