Traverse "curve" on a grid from a given cell

Question

I've got a binary grid which is mostly zero but contains a bunch of 1-cell wide lines (e.g. obtained from a ridge extraction). You can obtain an example input from the following line:

grid = ArrayPlot[
  Table[
    Sin[x + Pi/2*Tanh[y/3]] > -0.99 // Boole, 
    {x, -10, 10, 0.2}, 
    {y, -10, 10, 0.2}
  ], 
  Frame -> None, 
  ImageSize -> 101, 
  PlotRangePadding -> None
] // Binarize // Thinning // ImageData;

My actual data is a bit more complicated, but I've made sure that no white pixel has three white neighbours, so that the (Moore-)connected components really are just curves.

I've also got the data in the form of a list of coordinate pairs:

coords = Position[grid, 1];

Now I've got one of those points. It doesn't matter how I obtain it, but for testing you can just do:

startingPoint = RandomChoice[coords];

Now my goal is to walk along the line from that point in both directions either for a fixed number of steps or within a fixed (Euclidean or Chebyshev) distance (10 steps/cell widths, say) and collect all those points. That is, I want the (ordered) set of points which are on the same line as the startingPoint and in its immediate vicinity. I could either work with two ordered lists, one for the line segment on each side of startingPoint, or a single ordered list from one end of the line segment to the other.

I can easily solve this with some procedural code, where I just repeatedly look at the neighbours of the last point(s) I added and add their other neighbour to the list(s) as well. However, this doesn't seem like a very Mathematica-like solution.

What would be an efficient and idiomatic way to obtain this neighbourhood along the line?

@corey979 Unfortunately, the workings of FindCurvePath seem to be a bit too magical to use it reliably. It's not really clear when it tries to close a loop and when not, and when it suddenly breaks up the line. — Martin Ender, Mar 23 '17 at 10:07

george2079 · Answer 1 · 2017-03-23T19:56:22.973

4

a graph approach:

p = Position[grid, 1];
v = Range@Length@p;
e = UndirectedEdge @@@ (Select[Subsets[v, {2}], 
    EuclideanDistance @@ p[[#]] <= Sqrt[2] &]);
g = Graph[v, e, VertexCoordinates -> Reverse /@ p];
     HighlightGraph[g, Select[v, GraphDistance[g, 240, #] < 10 &], 
     VertexSize -> 2]

HighlightGraph[g, 
 Map[Style[#, Hue[RandomReal[{0, 2/3}]]] &, ConnectedComponents[g]], 
 VertexSize -> 2]

Edit: a faster way to generate the edge list.

e = Flatten@Outer[Function[{pos},
        If[pos != {},
         UndirectedEdge[#1, pos[[1, 1]]],
         Unevaluated@Sequence[]]]@
        Position[p, p[[#1]] + #2, 1] &, v,
        {{0, 1}, {1, -1}, {1, 0}, {1, 1}}, 1];

edited Mar 23 '17 at 19:56

answered Mar 22 '17 at 21:03

george2079

38,913
1
43
110

I like this because it lets me specify the distance as a number of steps and because it's easy to order the points along the line. – Martin Ender Mar 23 '17 at 10:14
I wonder whether it might be more efficient to generate the graph with RelationGraph instead of Subsets and Select. – Martin Ender Mar 23 '17 at 10:15
@MartinEnder Probably so. RelationGraph is "new in 10.2", I don't have it to try. NearestNeighborGraph might also be useful (also new in 10.2) – george2079 Mar 23 '17 at 14:13

Ali Hashmi · Answer 2 · 2017-03-24T02:02:12.863

With some credits to @C.E.

(* which trajectory the point lies on *)
whichTrajectory[grid_, startpos_] := With[{p = MorphologicalComponents[grid]},
With[{labels = (Union@Flatten@p)[[2 ;;]]},
{p, First @@Position[(Position[p, #]~MemberQ~startpos) & /@ labels, True]}]
];

(* to get nearest neighbours *)
ClearAll[createNearestFunc];
SetAttributes[createNearestFunc,HoldFirst];
createNearestFunc[f_][t:g_[x_,y_]] := f@Position[First@t, Last@t];

startingPoint = RandomChoice@coords; (* coords you already know *)

nf = createNearestFunc[Nearest] [whichTrajectory[grid, startingPoint]];
collections = nf[startingPoint, 10];

findOrder[pts_] := Module[{graph = NearestNeighborGraph[collections, 2],
subsets, euclideandist, periphery, path},

subsets = Subsets[GraphPeriphery[graph], {2}];
euclideandist = N@EuclideanDistance[##]&@@#&/@Subsets[GraphPeriphery[graph], {2}];
periphery = Part[subsets, First @@ Position[#, Max[#]] &@euclideandist];
path = MaximalBy[FindPath[graph, Sequence @@ periphery, Infinity, All], Length]
];


Show[ListPlot[coords, PlotStyle -> Black],ListPlot[collections, PlotStyle -> Red]]

New Edit: finding the ordered list of points on the trajectory. Note: the steps below have been packaged into the function findOrder above

We can find the ordered sequence of points on the trajectory using a graph.

suppose we have a collection of points

collections = {{79, 42}, {79, 41}, {78, 43}, {79, 40}, {78, 44},
{80, 39}, {77, 45}, {80, 38}, {77, 46}, {80, 37}}

then drawing a NearestNeighbourGraph for nearest 2 points

graph = NearestNeighborGraph[collections, 2, VertexLabels -> "Name"]

subsets = Subsets[GraphPeriphery[graph], {2}];

euclideandist = N@EuclideanDistance[##] & @@ # & /@Subsets[GraphPeriphery[graph], {2}];

peripherypos = First @@ Position[#, Max[#]] &@euclideandist;

periphery = subsets[[peripherypos]];

path = MaximalBy[FindPath[graph, Sequence @@ periphery, Infinity, All], Length]

(* {{{77, 46}, {77, 45}, {78, 44}, {78, 43}, {79, 42}, {79, 41}, {79,40},
{80, 39}, {80, 38}, {80, 37}}} *)

Now seeing the Line

Graphics[{Red, PointSize[Large], Point[#],Thick, Blue, Line[#]}] &@path

I just noticed that I haven't been particularly clear about the fact that I do need the points ordered along the line (either the two separate halves from the starting point, or the full line segment from one end to the other). This is a decent approach though. — Martin Ender, Mar 23 '17 at 10:10
@MartinEnder you can do SortBy[collections, Last] to sort them — Ali Hashmi, Mar 23 '17 at 11:37
@MartinEnder this will not work for horizontal case. Let me rethink about it ! — Ali Hashmi, Mar 23 '17 at 11:48
@MartinEnder please check the new edits. I guess this should work — Ali Hashmi, Mar 23 '17 at 23:22
@MartinEnder finally after all these edits i quit. i think its doing a good job :) — Ali Hashmi, Mar 24 '17 at 02:08

score 2 · Answer 3 · edited Apr 13 '17 at 12:56

Adapting this answer using ComponentMeasurements:

renumber = Module[{i = 1}, # /. 1 :> i++] &;
grid2 = renumber@grid[[;; -2, ;; -3]]; (* removed last row and last 2 columns*)

vertices = ComponentMeasurements[grid2, "Label"][[All, 1]];
centroids = ComponentMeasurements[grid2, "Centroid"];
neighbors = ComponentMeasurements[grid2, "Neighbors"];
edges = UndirectedEdge @@@ DeleteDuplicates[Sort /@ Flatten[Thread /@ neighbors]];

gr = Graph[vertices, edges, VertexCoordinates -> centroids]

Mathematica graphics

HighlightGraph[gr, {Style[150, Green], NeighborhoodGraph[gr, 150, 20]}, 
 VertexLabels -> {150 -> Placed["Name", Above]}, 
 VertexSize -> {_ -> .75, 150 -> 1.2}]

Mathematica graphics

Traverse "curve" on a grid from a given cell

3 Answers3