20
$\begingroup$

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.

Mathematica graphics

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?

$\endgroup$
  • $\begingroup$ ImageCorners reference.wolfram.com/language/ref/ImageCorners.html may be useful for step 1. $\endgroup$ – DavidC Apr 23 '15 at 11:52
  • $\begingroup$ For "more complicated" examples than the one presented or even complex networks ImageCorners results in way too many false positive findings. $\endgroup$ – Kardashev3 Apr 27 '15 at 10:41
16
+125
$\begingroup$

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]

enter image description here

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

enter image description here

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

enter image description here

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]

enter image description here

and segment that:

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ Doesn't MorphologicalTransform[img, "EndPoints"] help you with the first piece of code? $\endgroup$ – Dr. belisarius Jul 22 '15 at 18:22
  • $\begingroup$ @belisarius: yes, if I'd have known about it ;-) $\endgroup$ – Niki Estner Jul 22 '15 at 18:38
  • $\begingroup$ Your proposed method connects only endpoints. The idea should be to connect interrupted paths as well as generate new junctions. The images that have to be reconstructed need many of such junctions :(. $\endgroup$ – Kardashev3 Jul 27 '15 at 8:37
  • $\begingroup$ @Kardashev3: Then maybe you should have added a comment while the bounty was still active. Since you didn't, the obvious conclusion was you were either happy with the answer or lost interest. $\endgroup$ – Niki Estner Jul 27 '15 at 9:23
  • $\begingroup$ @nikie: I am very sorry about this. But I already mentioned the necessary forming of new branches in other comments. $\endgroup$ – Kardashev3 Jul 29 '15 at 15:36
4
$\begingroup$

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]

Mathematica graphics Mathematica graphics

$\endgroup$
  • $\begingroup$ Unfortunately this works for very small gaps only. Reconnection of larger interruptions or the forming of branches cannot be done using ListCurvePathPlot. $\endgroup$ – Kardashev3 Apr 27 '15 at 10:39

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.