I am not exactly working with blood vessels, but my problem is similar enough that I think someone who does work with blood vessels can answer this question for me.

I have several GPS points of a car driver going through a city. I would like to reconstruct the city map from these points.

The problem is that, unlike traditional vein extraction, there are intersections.

This is an example of the data I have: (reproduced as an image, but I have the actual coordinates.)

city map

This is ONE example of what I would like it to become. (I have only segmented some of the roads.) I would be happy however if the algorithm classified each intersection as the start of a new road.

city map2

I can and I have used image processing to reproduce a contiguous map (using dilation and skeletonization). The process is not perfect, but if the algorithm works only with lines, I can produce the map in that form.

Note: (1) the map is much bigger than this, and it is unfeasible to do it by hand. (2) And no, I cannot use a real map.

Related posts:

  • $\begingroup$ Can I please ask if there was any resolution to this question? $\endgroup$ – A_A Mar 12 '17 at 13:12
  • $\begingroup$ I have not answered yet because we are still experimenting with the data. I will reply very soon. But thank you very much for your invaluable contribution. $\endgroup$ – Ricardo Cruz Mar 13 '17 at 15:21
  • $\begingroup$ I have accepted your answer as the final solution. But please see my long answer below. Your answer was very helpful and inspirational, thanks. $\endgroup$ – Ricardo Cruz Mar 13 '17 at 15:41

Your problem is quite different than the vein extraction one. In fact, it is much mucg easier to solve with the "geometric" definition of the Hough Transform.

The Hough Transform is used exactly to detect lines in an image. It achieves this by integrating the brightness values of an image towards some angle $\theta$. So, for example, for angle $\theta = 0 ^ \circ$, the Hough Transform produces the sum of each row in the image, for angle $theta = 90 ^ \circ$, the Hough transform produces the sum of each column in the image...and so on.

The combined result of this is that points map to arcs, lines map to points and polygons map to specific configurations of points in the Hough Space.

The way to pick up a line in the image, in the Hough Space, is ultra simple. You just pick a maximum. But, that would give you an infinite line, not a line segment, which is what you have in your application.

The most straightforward way to pick up a line segment by the Hough Transform is to "mark" the pixels that are "responsible" for a particular accumulation. So, later on, when you pick up a maximum peak, you can "look" at the pixels that formed that sum and then connect them with a simple distance based rule. If distance is smaller than 4 pixels between any two pixels, then consider them as part of one line, otherwise create a new line.

For more information about all this, please see this link.

We now come to your case, you have two problems:

  1. Detect line segments from GPS points
  2. Connect line segments to form roads.

To solve #1, you can use the "geometric" definition of the Hough Transform WITHOUT having to work on the pixels of an image. Here is how this works:

  1. Find the centroid of all GPS points. Call this the centre of your region.
  2. Find the bounding box of all of your GPS points. These two steps define your "image space".
  3. Run a Hough Transform with the center and over the bounding box of the points. As you are shooting "rays" down each different bearing, make sure that you also store the points that are responsible for the sum. The sum is simple: Shoot a ray down a direction, whenever it hits a point increase a counter and store its location.
  4. Do a histogram of your Hough Space. It will help you in determining the threshold you have to apply so that you select the peaks that result in lines.
  5. Apply the Threshold by selecting peaks that are above the thrshold value. You now have a bunch of infinite strong lines in your hands.
  6. Look at the lists of points accumulated in each one of the directions and "integrate" them in lines if their distance is smaller than some threshold. You now have a bunch of line segments in your hands. These are your roads OR your road segments. The good thing about them is that they share a common point when they are connected (the junction). The bad thing about them is that they are unconnected when there is a "turn" or "curve".

You now have to solve problem number 2. Unfortunately, it is impossible to solve problem #2 without yet another application of a threshold. The philosophy here is simple:

  1. Pick up a line, call it $a$.
  2. Pick up all of the lines $a$ shares a common point with, call them $C$.
  3. Calculate all angles between $a$ and the lines in the set $C$.
  4. Merge two lines into one "road" if the angle is greater than some threshold. Obviously, if two lines form a $180^\circ$ angle, they are the same "road" but if they form a $90^\circ$ angle then they probably are not EXCEPT if set $C$ only has 1 element (1 line).

At the end of this process, you will have all connected segments that form roads.

It is not a difficult problem but it sure is a messy one.

Hope this helps.


My colleagues prefer that we do the work manually, which is probably the least effort solution.

The solution I was implementing is inspired by @A_A answer, but did not use any image processing technique.

Here it is in broad strokes:

For a point $i$,

  1. finds its neighbors within a radius $D$
  2. compute the normalize and the unormalized gradient of the neighbors relative to point $i$
  3. use the normalized gradient with the inverse of tangent to discretize in $A$ angles (for instance $A=4$ distinguishes the four directions: |, -, \ and /)
  4. for each of these angles, see if the point distribution is continuous. The idea is that points in a road will be very next each other and well distributed, while points that do not belong to the same road will have many blank spaces in the middle and are not as well distributed. (how we test for this, next)

Using this simple approach, we can tell without much error which points belongs to the same street as a given point $i$,


(the points in cyan were the ones chosen as being in the same road)

For some other point,


(the points in red were the ones chosen as being in the same road)

Point #4 was implemented in three different forms:

  • using Gini to test "inequality" of points
  • doing an 1D histogram (using a certain number of bins) of the gradient vectors and then testing it using a KS-test against the uniform distribution (the KS-test is defined as the maximum difference between the cumulative distributions)
  • the maximum difference between every two points in that angle.

This last test was based on the KS-test, and seems to be the test which has worked best.

The work is not complete however. My idea was now for each point to vote on what points belong to its road, and finish with a voting strategy approach. I was hoping this voting strategy, typical of ensembles, would fix some errors introduced by the previous technique. Such ensembles techniques are known to improve the robustness of underlying techniques.

But I did not find the time to finish it and, as mentioned, we are probably going with manual segmentation.

  • $\begingroup$ Thank you for your message. This is one of those exchanges i value a lot and keep me motivated to be coming back to this SE site. Just for the record, the proposed solution was not necessarily to be carried out over an image space, hence the term "Geometric". I have completely ignored the height element of the GPS which I suppose you are using here too. I have previously worked with "min dist" criteria but ended up with "webs" at the intersections. The Hough Trans modification was born out of a need to minimise that. $\endgroup$ – A_A Mar 13 '17 at 16:05

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.