Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a matrix of 0's and 1's forming a number of disjoint paths:


  RandomPaths
I would like to find the lengths of the paths, and from that "spectrum," the longest length (in the above example: 27, starting at {1,14}). (The shortest length possible is 3, just from how I generate these paths. There are never trees or cycles—just paths.) I can do it by identifying start 1-cells as either on the boundary, or having three 0-neighbors and one 1-neighbor, and then tracing from the start cell along the path. This is quite clunky.

Can anyone see a slicker method? Efficiency is not my main concern at the moment. Thanks for your ideas!

Here's the matrix displayed above:


{{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1,
   0, 0, 0, 0, 0, 0, 0}, {1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0,
   1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1}, {0, 1, 0, 0, 0, 0, 0,
   1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
  0, 1, 0, 0, 1, 0, 0, 1, 0}, {0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 
  1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1}, {0, 0, 0, 0, 1, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0,
   0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0,
   0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1}, {0, 1, 0, 0,
   1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 
  0, 0, 1, 0}, {0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
  0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0}, {1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 
  1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1}, {0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
   1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,
   0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, {1, 1, 0, 0, 1, 1, 0, 1, 0,
   0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1}, {0,
   0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 
  0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
  0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 1, 1, 1, 1, 0, 0, 
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 
  0}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
  0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
  1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0}, {0, 1, 1, 0, 1, 
  1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1,
   1, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
   0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
   0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 1,
   1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 
  1, 1, 1, 1}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 
  1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 
  1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0,
   1, 1, 1, 1, 1}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,
   0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0}, {1,
   1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 
  0, 1, 1, 0, 1, 1, 0}, {0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}
share|improve this question
    
You could use the image processing functions to find the connected components and sort them by length. The pixel positions can easily be mapped back to the matrix index –  rm -rf Jun 21 at 20:46
    
Possibly food for thought: blog.wolfram.com/2010/12/21/… and linked posts. –  Yves Klett Jun 21 at 20:57
1  
How about providing the matrix? –  Öskå Jun 21 at 21:00
    
@Öskå: Thanks for the suggestion; included now. –  Joseph O'Rourke Jun 21 at 21:29
1  
The answers attest to an amazing inventiveness of this community! –  Joseph O'Rourke Jun 22 at 21:07

8 Answers 8

up vote 25 down vote accepted

As I expressed in my comment above, it is possible (and easy) to use the image processing functions for this. Taking m to be the matrix above the following steps illustrate the idea:

img = Image@m;
ComponentMeasurements[img, "PerimeterCount"]
(* {1 -> 3, 2 -> 27, 3 -> 9, 4 -> 6, 5 -> 15, 6 -> 3, 7 -> 6, 8 -> 3, 9 -> 3, 10 -> 3, 
    11 -> 3, 12 -> 3, 13 -> 3, 14 -> 3, 15 -> 3, 16 -> 6, 17 -> 18, 18 -> 6, 19 -> 3, 
    20 -> 3, 21 -> 3, 22 -> 3, 23 -> 12, 24 -> 3, 25 -> 3, 26 -> 12, 27 -> 15, 28 -> 6, 
    29 -> 6, 30 -> 3, 31 -> 3, 32 -> 9, 33 -> 9, 34 -> 9, 35 -> 3, 36 -> 3, 37 -> 6, 
    38 -> 3, 39 -> 3, 40 -> 3, 41 -> 3, 42 -> 6, 43 -> 15, 44 -> 6, 45 -> 9, 46 -> 3, 
    47 -> 3, 48 -> 6, 49 -> 3} *)

From the above, the largest length is 27 (which you can also find programmatically) and the corresponding path is:

longestPath = SelectComponents[img, "PerimeterCount", # == 27 &]

To get the list of indices:

SparseArray[ImageData@longestPath]["NonzeroPositions"]
(* {{1, 14}, {2, 14}, {3, 14}, {4, 14}, {5, 14}, {5, 15}, {5, 16}, {5, 17}, {6, 17}, {7, 17}, 
    {8, 7}, {8, 8}, {8, 17}, {9, 8}, {9, 17}, {10, 8}, {10, 17}, {11, 8}, {11, 9}, {11, 10}, 
    {11, 11}, {11, 12}, {11, 13}, {11, 14}, {11, 15}, {11, 16}, {11, 17}} *)
share|improve this answer
    
Wow! Beautiful! –  Joseph O'Rourke Jun 21 at 21:30
    
To find the largest length you may use: Max[List @@@( ComponentMeasurements[img, "PerimeterCount"])[[;; , 2]]] –  Algohi Jun 22 at 8:08

Here is a Graph-based solution inspired by this Q&A where mat is your given matrix.

binaryGraph[mat_] := Module[{pos, edge, dedge}, pos = Position[mat, 1];
  edge = Select[Subsets[Range@Length@pos, {2}], 
    Last@# - First@# <= (Max@Dimensions@mat + 1) &];
  dedge = DeleteDuplicates[UndirectedEdge @@@ (Extract[edge, #] & /@ 
       With[{dist = N@(EuclideanDistance[pos[[#]], pos[[#2]]] & @@@ edge)}, 
         Flatten[Position[dist, #] & /@ DeleteDuplicates@N@Select[dist, # == 1 &]]])];
  Graph[dedge, VertexCoordinates -> Rule @@@ Thread[{Range@Length@pos, ({#2, -#1} & @@@ pos)}], 
    VertexLabels -> "Name", ImagePadding -> 20]]

g = binaryGraph@mat;
pos = Flatten@Position[#, Max@#] &@(VertexEccentricity[g, #] & /@ Range@Length@VertexList@g);
HighlightGraph[g, PathGraph[FindShortestPath[g, First@pos, Last@pos]],
  ImageSize -> 350]
Length@FindShortestPath[g, First@pos, Last@pos]

Mathematica graphics


If one wishes to go back on a matrix-like side, the exact path taken from the matrix mat can be given by posmax:

posmax = Position[mat, 1][[#]] & /@ FindShortestPath[g, First@pos, Last@pos];
newmat = mat;
(newmat[[Sequence @@ #]] = Red) & /@ posmax
ArrayPlot[newmat, Mesh -> All, MeshStyle -> Directive[{Thick, Gray}]]

Mathematica graphics

share|improve this answer
1  
I love that display of the matrix & paths--Very clever! –  Joseph O'Rourke Jun 21 at 21:59
    
@JosephO'Rourke Thank you :) I kinda like it too ;o) –  Öskå Jun 21 at 22:00
1  
@Öskå ++1 - esp. after your addendum –  eldo Jun 21 at 22:36
    
@eldo thanks :) the elements were here.., I just needed to put them together :P –  Öskå Jun 21 at 23:25

Assuming m is the test matrix; the following is a graph based solution:

nbours[x_, y_] := Module[{xr, yr, nts},
   xr = Select[{x - 1, x, x + 1}, # >= 1 && # <= Length@m &];
   yr = Select[{y - 1, y, y + 1}, # >= 1 && # <= Length@m &];
   nts = Cases[
     Table[{xrp, y}, {xrp, xr}]~Join~Table[{x, yrp}, {yrp, yr}], 
     Except[{x, y}]];
   If[m[[y]][[x]] == 1, {x, y} \[DirectedEdge] # & /@ 
     Select[nts, Part @@ {m}~Join~Reverse[#] == 1 &], {}]
   ];

map = Flatten[Table[nbours[x, y], {x, 1, Length@m}, {y, 1, Length@m}]];
g = Graph[map]; dist = GraphDistanceMatrix[g];
max = Max[Replace[dist, \[Infinity] -> -1, 2]];
VertexList[g][[#]] & /@ First@Position[dist, max]

Produces the pair of points that form the longest path of distance max.

share|improve this answer
1  
Nice use of various tools, e.g., GraphDistanceMatrix[]. –  Joseph O'Rourke Jun 22 at 1:46

Using a slightly modified version this answer in a related Q/A:

ClearAll[m, edges, v, vcs, g, cc, ap, ap2, pos];
m = Module[{i = 1}, mat /. 1 :> i++] (* where `mat` is the input matrix *);
v = ComponentMeasurements[m, "Label"][[All, 1]];
vcs = ComponentMeasurements[m, "Centroid"];
edges = UndirectedEdge @@@ DeleteDuplicates[
         Sort /@ Flatten[
           Thread /@ ComponentMeasurements[m, "Neighbors", CornerNeighbors -> False]]];

g = Graph[v, edges, VertexCoordinates -> vcs, ImageSize -> 400, 
             VertexSize -> .5, VertexStyle -> Blue, EdgeStyle -> Thickness[.01]];
cc = ConnectedComponents[g] 

Note: ConnectedComponents are sorted by length in version 9.0.1.0. For earlier versions you need to use SortBy[cc, -Length[#]&] (thanks: Oska)

ap = ArrayPlot[m, Mesh -> All, ImageSize -> 400, MeshStyle -> Directive[{Gray, Thick}],
      ColorRules -> Join[Thread[cc[[1]] -> Red], {0 -> White, _ -> Black}]];

ap2 = ArrayPlot[m, Mesh -> All, ImageSize-> 400, MeshStyle -> Directive[{Gray, Thick}],
      ColorRules ->  Flatten[{Thread[# -> 
       ColorData[{"DeepSeaColors", "Reverse"}][Length[#]/Max[Length /@ cc]]] & /@ cc,
            {0 -> White, _ -> Black}}, 2]] (* paths color-coded by length *);

Row[{HighlightGraph[g, cc[[1]], GraphHighlightStyle -> "Thick"], ap, ap2}]

enter image description here

To get the positions:

pos = First[Position[m, #]] & /@ # & /@ cc;
pos[[1]] (* positions of the longest-path elements *)
(* {{8, 7}, {8, 8}, {9, 8}, {10, 8}, {11, 8}, {11, 9}, {11, 10}, 
    {11, 11}, {11, 12}, {11, 13}, {11, 14}, {11, 15}, {11, 16}, 
    {11, 17}, {10, 17}, {9, 17}, {8, 17}, {7, 17}, {6, 17}, {5, 17},
    {5, 16}, {5, 15}, {5, 14}, {4, 14}, {3, 14}, {2, 14}, {1, 14}} *)
share|improve this answer
    
I was wondering how to tweak your code to make it work without the "diagonal neighbours", and it simply is CornerNeighbors -> False..! :D It's always obvious after seeing the answer :P –  Öskå Jun 22 at 16:16
    
@Öskå, it took me a while to remember CornerNeighbors:) –  kguler Jun 22 at 16:21
1  
I love the 3rd display using DeepSeaColors! –  Joseph O'Rourke Jun 22 at 16:24
    
@Öskå & Joseph, thank you both. Oska, I added a note on version differences. –  kguler Jun 22 at 16:34

Not nearly as clean as @rm -rf's solution, but here is another image processing solution. Once again m is taken to be the matrix above.

mc=MorphologicalComponents@m;

Max@Tally[Join @@ mc][[2 ;;, 2]]
(*27*)

If one wishes to vizualize the exact path:

p1=ArrayPlot[mc /. Rule @@@ Rest[Tally[Join@@mc]], Mesh -> All];
(* paths color-coded by length *)

p2=ArrayPlot[mc /. {Last@Commonest[Join@@mc, 2] -> Red, a_ /; a != 0 -> 1}, Mesh -> All];
(* colors longest path red *)

GraphicsRow[{p1, p2}]

enter image description here

share|improve this answer

Here is a direct approach w/o the graphics/graph functions.

Algorithm:

  • create an array with every '1' changed to a unique number (pad zeros around array to avoid dealing with edge issues )
  • loop over nonzero entries. If any of the four neighbors (on further thought you may only need to look left & down) swap indices to match - this is a global replacement over the whole matrix.
  • repeat the process until no unmatched neighbors are found (for this example it gets it on the first pass, but i think for some cases the iteration is needed )
  • the last step reassigns the region numbers so they are sequential and sorted by region size (not really needed depending on what you want to do with the result )

>

  (m = ArrayPad[m0 Partition[ Range[ Times @@ Dimensions[m0]] ,
   Last@Dimensions[m0]], 1];
  nonzero = Position[m, x_ /; x > 0, {2}];
  NestWhile[(Function[{p},
    If[ Length@# > 0 ,
        m = m /. m[[Sequence @@ p]] -> Min[#] ] &@ 
            Select[Table[ m[[Sequence @@ (p + s) ]] ,
               {s , {{0, 1}, {0, -1}, {-1, 0}, {1, 0}} }] ,
      (# != 0 && # != m[[Sequence @@ p]]) &] ] /@ nonzero; # ) & ,
      m , (#1 =!= #2 ) &, 2];
  m = m[[2 ;; -2, 2 ;; -2]] /. MapIndexed[ #1[[1]] -> First@#2 &,
        Reverse@
          SortBy[Tally@Cases[ Flatten@m , Except[0]], #[[2]] & ] ]) //Timing // First

0.0468

(surprisingly faster than ComponentMeasurements )

The result is a matrix with identifier number for each group (sorted by group size)

enter image description here

 nonzero = # - {1, 1} & /@ nonzero;(*shift because we removed the padding*)
 Graphics[{ 
    ColorData[60, "ColorList"][[Mod[m[[Sequence @@ #]] - 1, 21] + 1]],
     Disk[ {#[[2]], -#[[1]]}, .5]} & /@ nonzero]
     Reverse@SortBy[Tally@Cases[ Flatten@m , x_ /; x > 0], #[[2]] & ]

enter image description here

 Reverse@SortBy[Tally@Cases[ Flatten@m , x_ /; x > 0], #[[2]] & ]

{{1, 27}, {2, 18}, {5, 15}, {4, 15}, {3, 15}, {7, 12}, ...

share|improve this answer
    
I'm sure it has a name, but i've not formally studied the topic. At a glance It doesn't look exactly like the algorithms on that page, since I'm using mathematicas global pattern replace rather than iteratively flipping neighbors. –  george2079 Jun 23 at 20:03
    
OK, good job. +1. –  Pickett Jun 23 at 20:08
    
Impressive! I know this is supposed to be self-evident from the code, but I for one could use at least a high-level description of how this works... –  Joseph O'Rourke Jun 23 at 20:22
    
@george2079 - Great piece of code, but I don't like the lollipop-colors (It's difficult for me to immediately see the longest path). –  eldo Jun 23 at 22:24
    
Try this Hue[2/3 (m[[Sequence @@ #]]-1)/(Last@Union@Flatten@m]-1) instead of ColorData[] –  george2079 Jun 23 at 22:40

There are already several really good answers. I simply want to point out that this can be cast as using Graph in a way that is explicitly quite efficient. The steps are as follows.

(1) Make positions corresponding to 1 into vertices. (2) Find immediate neighbors and make all pairs into edges. (3) Create a Graph object with the above data and VertexCoordinates corresponding to the vertex coordinates (this is not quite the tautology it sounds like). (4) Find the connected components.

If we start with the original binary matrix then code to do this is as below. The reverse-transpose in the first step is to get the correct orientation for coordinates to correspond to the original picture.

smat = SparseArray[Map[Reverse, Transpose[mat]]];
verts = Most[ArrayRules[smat]][[All, 1]];

To get edges efficiently one might use a Nearest formulation. Alternatively one can use a loop approach that is effectively linear time in the number of vertices. With the slightly less efficient Nearest approach we will afterward remove duplicates where first vertex is lexicographically greater than second (since we'll be retaining the oppositely-ordered equivalents).

nf = Nearest[verts -> Thread[v[verts]]];
nbors = Map[nf[#, {5, 1}] &, verts];
edges = Flatten[
   Map[Thread[UndirectedEdge[First[#], Rest[#]]] &, nbors] /. 
    v -> Identity];
edges2 = Select[
   edges, #[[1, 1]] <= #[[2, 1]] && #[[1, 2]] <= #[[2, 2]] &];

Now we make this graph.

gr = Graph[edges];

Place the vertices where they belong. (Surely there must be a better way to do this in the graph creation step...)

PropertyValue[gr, VertexCoordinates] = VertexList[gr];

Extract the connected components.

paths = WeaklyConnectedComponents[gr];
pathlens = Map[Length, paths];
maxlen = Max[pathlens]

(* Out[279]= 27 *)

I will make no effort to compete with the far nicer pictures that have already been shown*.

HighlightGraph[gr, 
 Subgraph[gr, SelectFirst[paths, Length[#] == maxlen &]] ]

enter image description here

*Mostly because I'd lose. And badly.

share|improve this answer
    
Thanks for the annotations of your code! –  Joseph O'Rourke Jun 24 at 23:26

On the Wikipedia page for connected-component labeling a couple of algorithms are explained that probably are similar to what the component image processing functions are doing. Although I'm not really satisfied with the code, I'll go ahead and post an implementation of the two-pass algorithm, since it seems to be pretty popular for this task.

listNeighbors[matrix_, pos_] := DeleteCases[
  Extract[ArrayPad[matrix, 1], pos + {1, 1} + # & /@ {
     {1, 1}, {1, -1}, {1, 0}, {-1, 1}, {-1, -1}, {-1, 0}, {0, 
      1}, {0, -1}
     }], 0]

firstPass[matrix_] := Module[{linked, labels, nextLabel, neighbors, rows, cols},
  linked = {{}};
  labels = ConstantArray[0, {rows, cols} = Dimensions[matrix]];
  nextLabel = 1;

  Table[
   If[
    Extract[matrix, {i, j}] != 0,
    neighbors = listNeighbors[labels, {i, j}];
    If[
     Length[neighbors] == 0,
     linked = Insert[linked, nextLabel, {nextLabel, -1}];
     labels[[i, j]] = nextLabel;
     AppendTo[linked, {}];
     nextLabel++;,
     labels[[i, j]] = Min[neighbors];
     Do[
      linked[[label]] = Union[linked[[label]], neighbors],
      {label, neighbors}
      ]
     ]
    ],
   {i, rows}, {j, cols}
   ];
  {labels, linked}
  ]

secondPass[{labels_, linked_}] := Map[If[# == 0, 0, Min[linked[[#]]]] &, labels, {2}]

It can be used like this:

secondPass[firstPass[matrix]] // Colorize

components

The longest path can be retrieved by Tally as in some of the other answers.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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