8
$\begingroup$

I have this topological cadastral map, and others:

Cadastral map

It can be imported into Mathematica using

img = Import["http://i.stack.imgur.com/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.

$\endgroup$
  • 3
    $\begingroup$ Closely related: 19550 $\endgroup$ – C. E. Sep 17 '15 at 22:25
  • 1
    $\begingroup$ Did that link solve your problem? If so, please delete this post so that we don't have unanswered problems hanging around that actually already have solutions. If not, edit your post with more specific details about your problem. $\endgroup$ – march Sep 18 '15 at 3:26
  • 2
    $\begingroup$ possible duplicate of Image processing: Floor plan - detecting rooms' borders (area) and room names' texts $\endgroup$ – MarcoB Sep 18 '15 at 16:00
  • 1
    $\begingroup$ Not a duplicate, I aim to obtain a vector polygon, very different from the other answer $\endgroup$ – Zanfronio Sep 18 '15 at 16:04
  • 1
    $\begingroup$ Sidenote: I'm 99% certain that image processing is the wrong solution to this problem. The OpenStreetMap web interface has an "Export" button that exports a raw XML file (i.e. not a graphics file but a data file containing semantic data like "house" and "way" and longitudes/latitudes). Try to import that, everything else is a hack. $\endgroup$ – Niki Estner Sep 18 '15 at 22:07
9
$\begingroup$

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]]
        }
      }
   }]]

enter image description here

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]}]]

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ I understood about half of those words. Enough to know I can't help you with that... $\endgroup$ – Niki Estner Sep 19 '15 at 11:43
  • $\begingroup$ My english is really poor, sorry. I mean: the polygons created are a bit irregular. I'd like to fit their lines with those of the rectangle that better approximate the polygon's shape (if convex building), matching its lines to those of the road and the next building's. In other words, buildings should be aligned and squared as closer as possible to the original image. For non convex buildings, I guess two or more overlapping rectangles can be used. The aim is to get buildings aligned to each other and to the road, the rectangle's stuff is just an hypotetical tecnique. $\endgroup$ – Zanfronio Sep 19 '15 at 12:01
  • $\begingroup$ Sorry, I mistakenly deleted the previous comment, cannot recover. Its thanks are still actual! $\endgroup$ – Zanfronio Sep 19 '15 at 12:03
  • 1
    $\begingroup$ Sorry, that comment wasn't about your English (which is perfectly fine) but about my lack of knowledge regarding maps&projections. $\endgroup$ – Niki Estner Sep 19 '15 at 12:15
  • $\begingroup$ I think I can manage to project it someway. The main issue is what I outlined in the comment. Hope your knowledge in imageing-related issues may help $\endgroup$ – Zanfronio Sep 19 '15 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.