# Traverse "curve" on a grid from a given cell

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,
] // 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?

• Maybe FindCurvePath will be useful. Mar 22, 2017 at 17:33
• @corey979 That does look useful, thank you. Mar 22, 2017 at 17:35
• @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. Mar 23, 2017 at 10:07

a graph approach:

p = Position[grid, 1];
v = Range@Length@p;
e = UndirectedEdge @@@ (Select[Subsets[v, {2}],
EuclideanDistance @@ p[[#]] <= Sqrt &]);
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];

• 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. Mar 23, 2017 at 10:14
• I wonder whether it might be more efficient to generate the graph with RelationGraph instead of Subsets and Select. Mar 23, 2017 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) Mar 23, 2017 at 14:13

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] HighlightGraph[gr, {Style[150, Green], NeighborhoodGraph[gr, 150, 20]},
VertexLabels -> {150 -> Placed["Name", Above]},
VertexSize -> {_ -> .75, 150 -> 1.2}] 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. Mar 23, 2017 at 10:10
• @MartinEnder you can do SortBy[collections, Last] to sort them Mar 23, 2017 at 11:37
• @MartinEnder this will not work for horizontal case. Let me rethink about it ! Mar 23, 2017 at 11:48
• It's also problematic if the lines curve more strongly. Mar 23, 2017 at 12:26
• @MartinEnder please check the new edits. I guess this should work Mar 23, 2017 at 23:22