views:

422

answers:

4

I have a map that is cut up into a number of regions by borders (contours) like countries on a world map. Each region has a certain surface-cover class S (e.g. 0 for water, 0.03 for grass...). The borders are defined by:

  • what value of S is on either side of it (0.03 on one side, 0.0 on the other, in the example below)
  • how many points the border is made of (n=7 in example below), and
  • n coordinate pairs (x,y).

This is one example.

0.0300      0.0000           7
2660607.5   6332685.5   2660565.0   6332690.5   2660541.5   6332794.5 
2660621.7   6332860.5   2660673.8   6332770.5   2660669.0   6332709.5 
2660607.5   6332685.5

I want to make a raster map in which each pixel has the value of S corresponding to the region in which the center of the pixel falls.

Note that the borders represent sharp changes in S. The various values of S represent discrete classes (e.g. grass or water), and are not values that can be averaged (i.e. no wet grass!).

Also note that not all borders are closed loops like the example above. This is a bit like country borders: e.g. the US-Canada border isn't a closed loop, but rather a line joining up at each end with two other borders: the Canada-ocean and the US-ocean "borders". (Closed-loop borders do exist nevertheless!)

Can anyone point me to an algorithm that can do this? I don't want to reinvent the wheel!

+1  A: 

I'd recommend you to use a geometry algorithm library like CGAL. Especially the second example in the "2D Polygons" page of the reference manual should provide you what you need. You can define each "border" as a polygon and check if certain points are inside the polygons. So basically it would be something like

for every y in raster grid
  for every x in raster grid
    for each defined polygon p
      if point(x,y) is inside polygon p
        pixel[X][Y] = inside_color[p]

I'm not so sure about what to do with the outside_color because the outside regions will overlap, won't they? Anyway, looking at your example, every outside region could be water, so you just could do a final

    if pixel[X][Y] still undefined then pixel[X][Y] = water_value

(or as an alternative, set pixel[X][Y] to water_value before iterating through the polygon list)

schnaader
This would do the job, if all my borders were closed polygons... Any tips on how I can convert them to such?
JF
you cannot tell if you're inside or outside something which is not closed. you must describe your map using closed polygons, don't describe the borders between two items but the border of each item which will be a set of closed shapes.
fa.
Yes, I know, I realize that this is what I "must" do. I'm looking for existing algorithms that can do this, given the inputs that I have. See question for complete description of these inputs.
JF
Accepted this answer even though it doesn't quiiiite do the trick... Still the closest among available answers.
JF
A: 
  • first, convert all your borders into closed loops (possibly including the edges of your map), and indentify the inside colour. this has to be possible, otherwise you have an inconsistency in your data
  • use bresenham's algorithm to draw all the border lines on your map, in a single unused colour
    • store a list of all the "border pixels" as you do this
  • then for each border
    • triangulate it (delaunay)
    • iterate through the triangles till you find one whose centre is inside your border (point-in-polygon test)
    • floodfill your map at that point in the border's interior colour
  • once you have filled in all the interior regions, iterate through the list of border pixels, seeing which colour each one should be
Martin DeMello
A: 
  • choose two unused colors as markers "empty" and "border"
  • fill all area with "empty" color
  • draw all region borders by "border" color
  • iterate through points to find first one with "empty" color
  • determine which region it belongs to (google "point inside polygon", probably you will need to make your borders closed as Martin DeMello suggested)
  • perform flood-fill algorithm from this point with color of the region
  • go to next "empty" point (no need to restart search - just continue)
  • and so on till no "empty" points will remain
maxim1000
Doesn't this leave me with a finite (unfortunately non-zero) number of "border"-colored pixels?
JF
yes... it seems so. For each such point you can perform the same check as for points in 5-th step. Or you can even modify flood fill algorithm to check borders presented as analytical segments rather than sets of points
maxim1000
Heheh, I guess I "could". I just wish I could avoid reinventing the wheel. Surely someone out there somewhere has solved this task already.
JF
+2  A: 

The general case for processing this sort of geometry in vector form can be quite difficult, especially since nothing about the structure you describe requires the geometry to be consistent. However, since you just want to rasterize it, then treating the problem as a Voronoi diagram of line segments can be more robust.

Approximating the Voronoi diagram can be done graphically in OpenGL by drawing each line segment as a pair of quads making a tent shape. The z-buffer is used to make the closest quad take precedence, and thus color the pixel based on whichever line is closest. The difference here is that you will want to color the polygons based on which side of the line they are on, instead of which line they represent. A good paper discussing a similar algorithm is Hoff et al's Fast Computation of Generalized Voronoi Diagrams Using Graphics Hardware

The 3d geometry will look something like this sketch with 3 red/yellow segments and 1 blue/green segment:

sketch of 3d geometry

This procedure doesn't require you to convert anything into a closed loop, and doesn't require any fancy geometry libraries. Everything is handled by the z-buffer, and should be fast enough to run in real time on any modern graphics card. A refinement would be to use homogeneous coordinates to make the bases project to infinity.

I implemented this algorithm in a Python script at http://www.pasteall.org/9062/python. One interesting caveat is that using cones to cap the ends of the lines didn't work without distorting the shape of the cone, because the cones representing the end points of the segments were z-fighting. For the sample geometry you provided, the output looks like this:

voronoi rendering output

Theran
This is my second favorite answer. Thanks for taking the time to demonstrate!!
JF