Reconstruction of interrupted paths in images

Question

I want to reconstruct paths in images and have already implemented a working method in Mathematica. On the left side you can see an image with interrupted paths and on the right side the reconstructed paths.

In short, the method can be broken down into these steps:

1. Find all path endpoints and extrapolate the direction they are facing with splines (paths are assumed to be curved, not straight)
2. Put „search cones“ to the endpoints and let them all progress on their spline paths iteratively until:
3. A search cone meets another search cone (or hits an existing original or already reconstructed path), then stop the search and connect the two endpoints (or the endpoint with the point on the existing path) by interpolating a path using splines (this considers extrapolated endpoint directions and ensures more rounded paths)
4. Continue the iterative searching process until there are no more endpoints or a given search distance is exceeded

This methods allows reconstruction of paths and creation of branches (see images) at the same time. The code is parallelized where possible, sad thing is that the iterative searching process cannot be parallelized in its current implementation. The algorithm has to start over again after one interruption was reconstructed. This ensures that short interruptions are reconstructed first and longer interruptions (that might be reconstructed based on still interrupted paths) will be reconstructed in a plausible way.

The problem is that for large images with many interruptions the computation takes very long. Now I would like to know if there are other (maybe simpler) approaches to this problem. Or does Mathematica offer some methods that can be used to simplify the problem in parts or as a whole?

Answer 1

The result you want looks very close to a spanning tree - which can be found very fast. Maybe it needs some refinement, but I think it's a good first step.

Step 1: Find endpoints

img = Binarize@
  ImageTake[
   Import["http://i.stack.imgur.com/AdAT8.png"], {1, -1}, {1, 230}]

(* end points are points that are on one of the lines and have \
exactly 2 white pixels in their 3x3 neighborhood *)
endPointsImg = 
  Image@MapThread[
    Boole[#1 == 1 && #2 == 2] &, {ImageData[img], 
     ListConvolve[BoxMatrix[1], ArrayPad[ImageData[img], {1, 1}]]}, 2];

endPoints = PixelValuePositions[endPointsImg, 1];

HighlightImage[img, endPoints]

Using these, I can build a few useful data structures:

components = Image[MorphologicalComponents[img]];

componentIndices = Range[1., Max[ImageData[components]]];

componentIndexToEndPoints = 
  GroupBy[Thread[PixelValue[components, endPoints] -> endPoints], 
   First -> Last];

componentIndexToPoints = 
  AssociationMap[PixelValuePositions[components, #] &, 
   componentIndices];

componentIndexToNearest = Nearest /@ componentIndexToPoints;

Next, I'm going to approximate each curve using a bezier spline segment (Note I'm using my own control point -> bezier conversion, as BezierFunction doesn't seem to play well with FindMinimum):

(* bezierBasis . [controlpoints] gives points on the bezier spline *)
bezierBasis = 
  Transpose[
   Table[BernsteinBasis[3, k, Range[0., 1., 0.1]], {k, 0, 3}]];

approximateBezier[idx_] := 
 Module[{endPts, nearest, distanceToCurve, bezierControlPoints, 
   bezierPoints},
  endPts = componentIndexToEndPoints[idx];
  nearest = componentIndexToNearest[idx];

  distanceToCurve[pt_ /; VectorQ[pt, NumericQ]] := 
   Total[(pt - First[nearest[pt]])^2];

  bezierControlPoints = {First[endPts], {c1x, c1y}, {c2x, c2y}, 
    Last[endPts]};
  bezierPoints = bezierBasis.bezierControlPoints;

  bezierControlPoints /. Last@Quiet@FindMinimum[
      Total[
       distanceToCurve /@ bezierPoints], {{c1x, 
        bezierControlPoints[[1, 1]]}, {c1y, 
        bezierControlPoints[[1, 2]]}, {c2x, 
        bezierControlPoints[[-1, 1]]}, {c2y, 
        bezierControlPoints[[-1, 2]]}}]]

approximatedCurves = 
  AssociationMap[approximateBezier, componentIndices];
Show[img,
 Graphics[
  {Orange,
   {Thick, Dashed, BezierCurve /@ Values[approximatedCurves]},
   Table[Text[i, Mean[approximatedCurves[[i]]], {-1, -1}], {i, 
     componentIndices}]}]]

EDIT: Now that I think about it, approximating splines it probably way too complicated. You could get the curve directions at the endpoints from the 2nd order derivatives of the image just as easily, as I've done here

Finally, calculate "distances" between the end points and construct a spanning tree from that:

curveEnding[idx_, endPointIdx_] := 
 Module[{bezier = BezierFunction[approximatedCurves[idx]], 
   u = endPointIdx - 1.}, 
     {bezier[u], Normalize[D[bezier[x], x] /. x -> u]}]

curveDistance[{p1_, dir1_}, {p2_, dir2_}] := 
 SquaredEuclideanDistance[p1, 
   p2] + (dir1 . {{0, 1}, {-1, 0}} . dir2)^2*10

this distance is a sum of the squared euclidean distance between two points and the squared sine of the tangent angles (times some weight).

Now I construct a fully connected undirected graph with the end points as vertices, where the end points of one segment are connected with weight 0 (so they're guaranteed to be in the spanning tree) and the other end points are weighted with curveDistance

endPointDistances = 
  Table[Property[
    UndirectedEdge[node[{componentIndex1, endPointIndex1}], 
     node[{componentIndex2, endPointIndex2}]], 
         EdgeWeight -> 
     curveDistance[curveEnding[componentIndex1, endPointIndex1], 
      curveEnding[componentIndex2, endPointIndex2]]], 
       {componentIndex1, componentIndices}, {componentIndex2, 
    Select[componentIndices, #1 > 
       componentIndex1 & ]}, {endPointIndex1, 1, 2}, 
       {endPointIndex2, 1, 2}]; 

connectEndPoints = 
  Table[Property[
    UndirectedEdge[node[{componentIndex, 1}], 
     node[{componentIndex, 2}]], EdgeWeight -> 0], 
       {componentIndex, componentIndices}]; 

spanningTreeEdges = 
 EdgeList[FindSpanningTree[
   Graph[Flatten[{endPointDistances, connectEndPoints}]]]];

Show[img, 
 Graphics[{Red, 
   spanningTreeEdges /. {UndirectedEdge[node[{c1_, e1_}], 
       node[{c2_, e2_}]] :> 

      If[c1 != c2, 
       Line[{componentIndexToEndPoints[c1][[e1]], 
         componentIndexToEndPoints[c2][[e2]]}]]}}]]

Like I said, not exactly what you wanted, but close.

One last note: I'm not sure how large you images are. The graph is fully connected, so it has O(n^2) nodes. One way to reduce that number is to connect only "neighboring" curves.

For that, you can construct a distance transform of you your image:

HighlightImage[ImageAdjust[DistanceTransform[ColorNegate@img]], img]

and segment that:

watersheds = 
  WatershedComponents[DistanceTransform[ColorNegate@img], 
   Method -> "Basins"];
HighlightImage[watersheds // Colorize, img]

From there, you can easily find out which curve segments are close enough that they might be connected:

ComponentMeasurements[watersheds, "Neighbors"]

{1 -> {2, 3}, 2 -> {1, 3, 4, 6}, 3 -> {1, 2, 6}, 4 -> {2, 5, 6}, 5 -> {4, 6, 7}, 6 -> {2, 3, 4, 5, 7}, 7 -> {5, 6}}

But for the 7 segments in your test image, that would be complete overkill, so I didn't bother.

Answer 2

ListCurvePathPlot[ ] might work where the directions and distances are within its limits:

i = Thinning[Binarize@Import@"http://i.stack.imgur.com/CKwCQ.png"]
pts = Position[ImageData[i], 1, {2}];
Framed@ListCurvePathPlot[pts, Axes -> False]