6
$\begingroup$

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;

enter image description here

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?

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

3 Answers 3

4
$\begingroup$

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]

enter image description here

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

enter image description here

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];
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Mar 23, 2017 at 10:14
  • $\begingroup$ I wonder whether it might be more efficient to generate the graph with RelationGraph instead of Subsets and Select. $\endgroup$ Mar 23, 2017 at 10:15
  • $\begingroup$ @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) $\endgroup$
    – george2079
    Mar 23, 2017 at 14:13
2
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Mar 23, 2017 at 10:10
  • $\begingroup$ @MartinEnder you can do SortBy[collections, Last] to sort them $\endgroup$
    – Ali Hashmi
    Mar 23, 2017 at 11:37
  • $\begingroup$ @MartinEnder this will not work for horizontal case. Let me rethink about it ! $\endgroup$
    – Ali Hashmi
    Mar 23, 2017 at 11:48
  • $\begingroup$ It's also problematic if the lines curve more strongly. $\endgroup$ Mar 23, 2017 at 12:26
  • $\begingroup$ @MartinEnder please check the new edits. I guess this should work $\endgroup$
    – Ali Hashmi
    Mar 23, 2017 at 23:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.