4
$\begingroup$

I am trying to create an isochrone based on the OpenStreetMap data set. Everything works fine, I extracted data, processed it into a DAG, ran a Dijkstra algorithm over it. The result is a subset of the complete graph. This is an impression of the covered parts of the subset displayed over Google Maps:

all vertices of graph

However, when the area gets larger, the number of reached vertices gets very large and displaying like this gets slow. What I would like to do is convert the set of edges and vertices into a polygon. Basically, this should be posible by removing all of the inner edges, leaving just the edges around the boundary of the area and the edges pointing out from it. I know coordinates for all vertices and approximating each edge as a line would be fine. Larger inner areas should become holes inside the polygon.

My first attempt was to use an geospatial library (in my case the SqlServer spatial extensions), create a multiline from all of the edges and doing an ST_Buffer on it. Turns out to be very slow and memory consuming for large numbers of edges (> 1000)

I was thinking along the lines of finding small polygons in the set (turning left at every turn?) and removing every edge that is part of 2 of these polygons.

Extra image to use in comment below: sample graph

| cite | improve this question | |
$\endgroup$
  • $\begingroup$ Would convex hulls help? $\endgroup$ – AJed Dec 18 '12 at 19:06
  • $\begingroup$ @AJed: no, a convex hull will not do. An isochrone typically has more of a star sharp. You would need an alpha shape, but that is a very intensive computation. $\endgroup$ – Teun D Dec 18 '12 at 22:37
  • 1
    $\begingroup$ Thanks all. I received three responses, all of which where interesting and valuable. I accepted Marcs answer mostly because I learned the most from it. $\endgroup$ – Teun D Jan 8 '13 at 17:43
5
$\begingroup$

I think you can do this with marching squares, which is both fast and easy to implement.

To use marching squares you need to construct a scalar grid over your domain. Impose a regular grid over your region of interest. For each point in the grid, compute the distance (or time, etc.) between the grid point and your point of interest using a shortest path algorithm. Store this value at that point to serve as the scalar values of your grid. Then apply the marching squares algorithm to generate an isocontour for some isovalue $d$; all points on this contour will be distance $d$ away from your point of interest. More on the details of this algorithm on Wikipedia.

Additionally you can speed up the shortest path computation since your graph is a DAG. First topologically sort your DAG and then relax each edge, like you would in Bellman-Ford, in the order given by the sorting. The correctness follows from relaxing each edge along every path from the source node to every other node. Thus the single source shortest path can be computed for a DAG in $O(V+E)$.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks a lot for your answer. This approach is actually very much like something I have already tried. It works, but to get a satisfying resolution, the squares had to be quite small and the resulting shape still becomes unpractical. It also feels inefficient, as I already have all of the vectors of a boundary in my set. Pixelating feels like an unnecessary step. But it is the fastest solution I have been able to create so far. $\endgroup$ – Teun D Dec 19 '12 at 6:49
  • $\begingroup$ On a second look at the wikipedia page you referred, I must say that this algorithm does generate a nicer shape around the squares than I did. Must read further. $\endgroup$ – Teun D Dec 19 '12 at 6:55
4
$\begingroup$

For each vertex $i$, note its $(x_i,y_i)$ location and the distance $d_i$ from the source. Building the Voronoi diagram or the Delaunay triangulation of the pointset can be useful (this would not use $d_i$ data).

If you want something simpler and self-contained, then for each point find its closest neighbors (this info can also be read from the Delaunay triangulation) and compare their $d_i$ values. In particular, if the current point $i$ has no neighbors $j$ with $d_j>d_i$, then $i$ is on the boundary. Otherwise, consider each closest neighbor $j$ with $d_j > d_i$ and make decisions based on the straight-line (or street) distance $d_{ij}$. For example, $i$ would still be on the boundary if for all neighbors $j$ you find $d_i + d_{ij} > d_j$.

Depending on the exact desired output, you may want to add some additional neighboring vertices to the boundary identified by the above condition.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ And with $d_i$ you refer to the distance over the graph, not the distance in the Euclidian space, I guess? $\endgroup$ – Teun D Jan 2 '13 at 13:27
  • $\begingroup$ And once I have determined all of the boundary edges, do you have any pointers to how to make a polygon with potentially holes and islands out of this set? $\endgroup$ – Teun D Jan 2 '13 at 13:33
  • $\begingroup$ I was originally thinking about the Euclidean distance, but perhaps, you can use either. Once you determined the boundary edges/vertices, consider vectors emanating from the center to the boundary vertices. For each vector indexed $j$, find its angular coordinate $\rho_j$ and sort all vectors by $\rho_j$. This ordering should give you a star-like (not necessarily convex) polygon. I imagine you may need to perform some adjustments for your application, but hope the general idea is what you needed. $\endgroup$ – Igor Markov Jan 2 '13 at 18:26
  • $\begingroup$ One more point. Check the image of a sample graph I have added to the end of the question. I think your approach would wrongly identify vertex A as being not on the boundar (it is connected to another point with greater $d$), while it would wrongly say that B is on the boundary. Using d as distance over the graph that is. Using Euclidian distance, A and C1, C2 would still by excluded from the boundary. Or am I still not quite understanding your point? $\endgroup$ – Teun D Jan 9 '13 at 12:22
  • $\begingroup$ C1 and C2 can be added by special-case processing. As for A, consider this -- connect the vertex on the mid-left to the vertex on the lower right. This would seem to remove A from the boundary. So, whether A is on the boundary or not, cannot be decided by looking only at its neighbors. One approach to consider is to first find parts of the boundary using the algorithm I suggested, and then incrementally fill in the missing pieces. For example, once you find "extreme points" of the boundary, you could try connecting them with paths that traverse boundaries. $\endgroup$ – Igor Markov Jan 10 '13 at 22:50
3
$\begingroup$

I would suggest the following. The data you extracted from OSM is not enough. A DAG only represents the combinatorial data, but not the topological data. What you should do is to sort for every vertex of the DAG the emanating edges in cyclic order, say ccw. Also I would forget about the directions of the edges. With this information you can start at any boundary edge and then you simply walk along the outer face until you reach your starting position.

Since your graph is planar, the sorting takes just $O(n \log n)$.

It also seems that yo have vertices of degree 1 on the boundary. You should get rid of them, for example by doubling the incident edges.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I have some properties attached to both the vertices (lat, lng) and the edges (length, maxspeed, oneway direction). So I do have a bit more info than just the DAG and I actually have all of the information needed for your approach. How would you handle holes inside the area? And in your O(nlogn) estimate, is n the avarage number of edges from a vertex? Wouldn't the full procedure be O(square root of vertices), assuming that the number of contained vertices grows with ^2 compared to the number in the boundary? $\endgroup$ – Teun D Dec 19 '12 at 14:13
  • $\begingroup$ Since your graph is planar you have less than 3 times as many as edges than vertices. So the sum of the vertex degrees is $O(n)$, for $n$ the total number of vertices. Since $n\log n$ is a convex function, the worst case running time for sorting everything is $O(n\log n)$. $\endgroup$ – A.Schulz Dec 19 '12 at 14:47

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.