Detect buildings in cadastral map as polygons

Question

I have this topological cadastral map, and others:

It can be imported into Mathematica using

img = Import["https://i.sstatic.net/3udXJ.jpg"]

I want to retrieve a description of the buildings, to be uploaded to OpenStreetMap.

I'd like to extract features from it (edges, etc.) and get the buildings' shapes in the form of polygons. I tried EdgeDetect - it does get edges and borders, but outputs a raster image. What I want is a vector image instead, with lines for each edge of the building.

None of the functions I used worked for me (such as Export to SVG), and I've run out of ideas.

What should I try? How do I best fit polygons to the result of EdgeDetect?

19550 is a (not very) similar question with answers that may be useful.

Answer 1

If you know that all the areas are convex, then this is quite easy: ComponentMeasurement has a measurement "ConvexVertices" that returns just that: the vertices of the convex hull:

bin = Binarize[img];
comp = ComponentMeasurements[
   DeleteSmallComponents[bin], {"Area", "ConvexVertices", 
    "ConvexCoverage"}, 50 < #1 < 10000 &, CornerNeighbors -> False];

the vertices are unordered (or maybe in the order they were discovered?), so I need to sort them:

sortClockwise = 
  Function[pts, 
   With[{c = Mean[pts]}, SortBy[pts, ArcTan @@ (# - c) &]]];

Then I can display them:

colors = ColorData[97];
Show[bin, Graphics[
  {
   EdgeForm[Red],
   comp /. 
      {(idx_ -> {area_, pts_, convexCover_}) :> 
        {If[convexCover > .75, 
            Directive[colors[idx], Opacity[0.5]], 
            Red],
        Polygon[sortClockwise[pts]]
        }
      }
   }]]

I've marked the polygons with low "ConvexCoverage" red: those are the non-convex buildings. If you can live with the convex hull for these, you're done.

If you absolutely need the vertices for concave buildings, things get a little uglier. Getting the border vertices for each component is easy enough: You can get a Mask for each component, that's a binary mask for the component:

compConvex = 
  ComponentMeasurements[
   DeleteSmallComponents[bin], {"Area", "ConvexCoverage", "Mask"}, 
   50 < #1 < 10000 && #2 < 0.75 &, CornerNeighbors -> False];

The idea is then to use MorphologicalPerimeter to get the perimeter of said mask, and use PixelValuePositions to convert the resulting binary image to pixel coordinates. Sadly, that gives us a very "unclean" perimeter, because of the label texts in you image. I've tried to remove those using a Closing morphological filter.

Sorting the vertices of a convex polygon is more difficult, too. I'll use FindShortestTour as a quick&dirty way to sort them (that's an approximative algorithm to an NP complete problem, so you're not guaranteed a perfect solution, though!)

outlineFromComponent[idx_ -> {area_, convexCover_, mask_}] := 
  Module[{vertices},
   vertices = 
    PixelValuePositions[
     MorphologicalPerimeter[Closing[Image[mask], DiskMatrix[5]]], 1];
   vertices = vertices[[FindShortestTour[vertices][[2]]]]];

Now we get the "right" outlines for the convex areas, too.

Show[img, 
 Graphics[{EdgeForm[Red], Opacity[0.2], Red, 
   Polygon[outlineFromComponent /@ compConvex]}]]

Note: The polygons still contain one vertex for each boundary pixel. You might have to implement something like the Ramer-Douglas-Peucker algorithm to reduce the number of vertices.