Finding connected components in an array

Question

I have an $w\times h$ array with non-negative integers.

I wish to partition the set of coordinates $(c,r)$ with $1\leq c \leq w$ and $1\leq r \leq h$, in subsets, so that each subset corresponds to a cluster of the same value in the array.

Two coordinates are neighbouring if they either are on top of another, or one is down-left of the other. This corresponds to finding clusters in a parallelogram.

belisarius was nice to provide a picture:

I have some code to find one connected cluster in an array (I call it gt):

FloodFillExtractTile[gt_, {sr_, sc_}] := 
 Module[{r, c, toExplore, visited = {}},
  toExplore = {{sr, sc}};
  While[Length@toExplore > 0,
   (* Pop *)
   {r, c} = Last@toExplore;
   toExplore = Most[toExplore];
   AppendTo[visited, {r, c}];
   (* Down-left *)

   If[c > 1 && r < h && 
     gt[[r + 1, c - 1]] == gt[[r, c]] && ! 
      MemberQ[visited, {r + 1, c - 1}],
    AppendTo[toExplore, {r + 1, c - 1}];
    ];
   (* Down-right *)

   If[ r < h && 
     gt[[r + 1, c]] == gt[[r, c]] && ! MemberQ[visited, {r + 1, c}],
    AppendTo[toExplore, {r + 1, c}];
    ];
   (* Up-right *)

   If[r > 1 && c < w && 
     gt[[r - 1, c + 1]] == gt[[r, c]] && ! 
      MemberQ[visited, {r - 1, c + 1}],
    AppendTo[toExplore, {r - 1, c + 1}];
    ];
   (* Up-left *)

   If[ r > 1 && 
     gt[[r - 1, c]] == gt[[r, c]] && ! MemberQ[visited, {r - 1, c}],
    AppendTo[toExplore, {r - 1, c}];
    ];
   ];
  Return@visited;
];

Then FloodFillExtractTile[gtp, {1, 1}] gives the component connected to the upper left hand corner. However, this method feels ugly, and extending it to all components feels even uglier.

I was looking at Gather, but the problem is that it wants all points in a cluster to be equal, see for example Gather dependency on list order?

Edit: So this is the type of arrays I am looking at. The $6$, the $3$'s, the $2$'s and the $0$'s are in one component respectively, but there are two components with $1$'s.

Being neighbours means being adjacent down-left, down-right, up-left and up-right.

Now, the rows is stored just like regular rows in a rectangular matrix, so that is why this translates to a bit strange criteria for being neighbours.

This is the best code I have so far, first extract points with the same value, then do connected-component analysis on those parts.

GetGTTiles[gtp_] := 
  Module[{testSame, testEdge, h, w, pts, sameClusters, getEdges, 
    tiles},
   {h, w} = Dimensions[gtp];
   pts = Join @@ Table[{r, c}, {r, h}, {c, w}];
   testSame[{r1_, c1_}, {r2_, c2_}] := (gtp[[r1, c1]] == 
      gtp[[r2, c2]]);
   testEdge[{r1_, c1_}, {r2_, 
      c2_}] := (gtp[[r1, c1]] == 
       gtp[[r2, c2]]) &&
     ((c1 == c2 && 
         Abs[r1 - r2] <= 1) || (c1 == c2 - 1 && 
         r1 == r2 + 1) || (c1 == c2 + 1 && r1 == r2 - 1));
   sameClusters = Gather[pts, testSame];
   getEdges[clust_] := 
    Join @@ Outer[If[testEdge[#1, #2], #1 -> #2, Sequence @@ {}] &, 
      clust, clust, 1];
   tiles = 
    Join @@ (ConnectedComponents[Graph@getEdges[#]] & /@ sameClusters);
   Return@tiles;
   ];

This is the output I expect for the example given (5 clusters found):

{{{1, 1}}, {{1, 2}, {2, 1}}, 
{{1, 3}, {2, 3}, {3, 2}, {4, 1}, {3, 1}, {2, 2}}, 
{{1, 4}}, 
{{2, 4}, {3, 4}, {4, 4}, {5, 4}, {5, 3}, {4,3}}, 
{{3, 3}, {4, 2}, {5, 2}, {5, 1}}}

EDIT: So this is the final code, based on belisarius solution:

GTTiles[gtp_List] := Module[{fromEuclidean, toEuclidean,
    getOneTile, elements, elmPos, pts, tile, tiles},

   (* This is used to changefrom different coordinate systems. *)

   fromEuclidean[{r_, c_}] := {r, (c - r)/2 + 1};
   toEuclidean[{r_, c_}] := {r, 2 c + r - 2};

   getOneTile[pts_List, maxDist_?NumericQ] := Module[{f},
     f = Nearest[pts];
     FixedPoint[
      Union@Flatten[f[#, {Infinity, maxDist}] & /@ #, 1] &, {First@
        pts}]];

   elements = Union @@ gtp;
   elmPos = (toEuclidean /@ Position[gtp, #]) & /@ elements;
   (* This is really strange code. *)
   tiles = Flatten[Flatten[
        Reap[NestWhile[Complement[#,
             Sow@getOneTile[#, N@Sqrt@2]] &, #, # != {} &]][[2]], 
        1] & /@ elmPos, 1];
   tiles = Map[fromEuclidean, tiles, {2}];
   Return@tiles;
   ];

Answer 1

l1 = {{6, 3, 2, 1}, {3, 2, 2, 0}, {2, 2, 1, 0}, {2, 1, 0, 0}, {1, 1, 0, 0}};
l = Riffle[#, ""] & /@ l1;
els = Union @@ l;
par = MapIndexed[PadRight[PadLeft[#1, #2[[1]] + Length@#1 - 1, ""], 
                          Length@l + Length@#1 - 1, ""] &, l];
eachElm = Position[par, #] & /@ els;
getOneCluster[pts_List, maxDist_?NumericQ] :=
  Module[{f},
   f = Nearest[pts];
   FixedPoint[Union@Flatten[f[#, {Infinity, maxDist}] & /@ #, 1] &, {First@pts}]];
clusters = 
  Flatten[Flatten[
      Reap[NestWhile[
         Complement[#, 
           Sow@getOneCluster[#, N@Sqrt@2]] &, #, # != {} &]][[2]], 1] & /@ eachElm, 1];
Grid[par, 
 ItemStyle -> {Automatic, Automatic, Flatten@MapIndexed[#1 -> Hue[#2[[1]]/3] &, 
                                                       clusters, {2}]}]

Answer 2

I've added this as a separate answer, since the method is so different from my other. I pondered this interesting problem today and came up with the following (Edit - see latter part of post for a much faster alternative...):

weirdNeighbors3[array_, output_: "components"] :=

 Module[{t = 0, cnt = 1, dims = Dimensions[array], 
   lt = Rest[array][[All, ;; -2]], rt = Most[array][[All, 2 ;;]], 
   check, a, ruls, clusts},

  clusts = 
   Transpose[Flatten /@ MapAt[Replace[#, _ -> cnt++, 1] &, 
      Split /@ Transpose[array], {All, All}]];

  check = Boole[MapThread[SameQ, {lt, rt}, 2]];

  While[t =!= {},

   t = DeleteCases[
     Transpose[{Flatten[Pick[Rest[clusts][[All, ;; -2]], check, 1]],
       Flatten[Pick[Most[clusts][[All, 2 ;;]], check, 1]]}], {} | {a_,
        a_}, 1];

   If[(t = 
       DeleteDuplicates[Flatten[ReplacePart[#, {a_, 1} /; 
              a > 1 && a <= Length[#] :> #[[a - 1, 2]]] & /@ 
          GatherBy[t, First], 1]]) === {}, Continue[]];

   clusts = 
    clusts //. Dispatch[
      Rule @@@ (Module[{f, g}, g[_] = True; f[{x_, _}] := (g[x] = False; True);
         Cases[t, {_, _?g}?f]])];

   ];
  Switch[output, "tuples", 
   GatherBy[Flatten[Array[{##} &, dims], 1], Extract[clusts, #] &], 
   "raw", clusts, "components", ArrayComponents[clusts]]
  ]

A brief explanation: The array is sliced by column, and initial cluster "IDs" are assigned, with any consecutive equivalent column members getting the same ID. A selection array is then built to allow fetching of diagonal pairs for checking of diagonal chains. Finally, in the While, rules are generated and massaged based on what is currently in the "clusts" cluster tracking array, and are then "percolated" through the array. This continues until there are no updates to be made (FixedPoint could be used here just as well).

The way the initial slicing is done (direction and order, e.g.) and the way the rules are created and "massaged" can be trivially changed to reflect other concepts of "connected", or to enhance performance for each case (e.g., there's another 20% in performance by doing the OP case by diagonals, I felt the code easier for readers to grok by going linear).

The default output is as MMA's ArrayComponents which I thought most useful to readers needing a "personalized" ArrayComponents. The optional argument using "tuples" outputs the array coordinate chains as per the OP (n.b.: the sort order differs from the OP example, due to the technique, it is of course trivial to re-sort). There is also a "raw" output, this is mainly for my debugging but gives the results as "massaged" by the process.

Some timing comparisons. The "normal" test used arrays generated with RandomInteger[{0, 3}, {dim, dim}], providing a reasonably "dense" test case with good probability of varying chain sizes. The "down" test used

makearrayD[dims_] := Module[{toprow = RandomInteger[dims[[1]]*5, dims[[2]]]},
  Table[toprow - 2 row, {row, 1, dims[[1]]}]]

to produce the "row-decreasing" arrays of the OP. Timings shown are relative performance, same with the extreme removed so details can be seen, overall performance trend in Log format, and performance advantage. As usual, timings on the cigar-lounge netbook...big thanks to belisarius and Mr. Wizard for enlightening ideas that got an extra 10% boost.

Update: Well, the update that I've now replaced below that was quite a bit faster than my other soultions and orders of magnitude faster than the AA on complex arrays just got blown out of the water in the same way.

The following was so fast even on my cigar-lounge netbook, I thought it broken and checked results several times before I believed it. I think this is about as much blood as can be squeezed from this turnip:

weirdNeighbors5[arg_] := Module[{dims = Dimensions[arg], cd, dd, p},

  cd = {#, # + {1, 0}} & /@ SparseArray[Unitize[Most[arg] - Rest[arg]], Automatic, 1][
     "NonzeroPositions"];

  dd = {# + {0, 1}, # + {1, 0}} & /@ 
    SparseArray[Unitize[Most[arg][[All, 2 ;;]] - Rest[arg][[All, ;; -2]]], 
      Automatic, 1]["NonzeroPositions"];

  Switch[Length[p = Join[cd, dd]],
      0, {{}, Flatten[Array[{##} &, dims], 1]}, 
      Times @@ dims, {p, {}}, _,
   {ConnectedComponents@Graph[UndirectedEdge @@@ p], 
    Complement[Flatten[Array[{##} &, dims], 1], Flatten[p, 1]]}]]

Performance snippet:

In cases where there are no connected components, the advantage is blinding: on a 150X150 array, timings were 313 seconds vs 0.0156 seconds - an over twenty-thousand to one advantage. Nice when there's a need to just determine if there are in fact any connected components (where I'm using it with a different rule for "neighbors").

Answer 3

Might be useful even though there's an accepted answer - about five to twenty-five times faster in my minor tests:

ClearAll[neighs, trace, weirdNeighbors]

(* get possible "neighbors" helper function *)
neighs[arr_, place_] := Module[{dims = Dimensions[arr], n},
   n = DeleteCases[place + # & /@ {{1, 0}, {-1, 0}, {1, -1}, {-1, 1}},
     {a_, b_} /; a == 0 || b == 0 || a > dims[[1]] || b > dims[[2]]]];

(* trace potential "neighbor" paths helper funtion *)
trace[list_, ele_] := Module[{results = {ele}},
   If[MemberQ[list, ele + {1, 0}], results = Join[results, trace[list, ele + {1, 0}]]];
   If[MemberQ[list, ele + {1, -1}], results = Join[results, trace[list, ele + {1, -1}]]];
   DeleteDuplicates[results]];

(* do the work function *)
weirdNeighbors[array_] := 
 Module[{local, td = Dimensions[array], ta, localc, prelim, tu, got, ss, fu, sets, f},

  local = array /. (0 -> Max[array] + 1);
  ta = Flatten[Array[{##} &, td], 1];
  localc = ArrayComponents[local];
  prelim = GatherBy[ta, localc[[Sequence @@ #]] &];
  tu = Union @@ localc;

  got = {ss = Select[prelim[[#]], 
        Function[arg, MemberQ[Extract[localc, neighs[localc, arg]], 
          Extract[localc, arg]]]], Complement[prelim[[#]], ss]} & /@ tu;

  fu = Function[arg, trace[got[[arg, 1]], #] & /@ got[[arg, 1]]] /@ tu;

  sets = Function[arg, Union @@@ Gather[arg, Intersection[#1, #2] != {} &]] /@ fu;

  Transpose[{tu, Flatten[array] // DeleteDuplicates, got[[All, 2]], sets}]]

Output is lists each consisting of the component number, what element that component corresponds to, a list of isolated members of that component, and a list of connected sets of that component. Results are consistent with my reading of the slightly ambiguous OP:

test = RandomInteger[{0, 5}, {5, 5}];

{time, result} = Timing[weirdNeighbors[test]];

Column[{time, test // MatrixForm, result}, Left, 2]

I've not done super extensive testing/timing, but it appears to do pretty well. Using a test of RandomInteger[{0, 20}, {100, 100}], it completed in 6.69 seconds. The same test with the accepted answer took 173.13 seconds, or about twenty-five times slower. With a test array of RandomInteger[{0, 100}, {200, 200}], it finished in 16.24 seconds, the AA took over 25 minutes, about ninety-five times longer. Caveat - as usual, timings on my cigar-lounging netbook, expect order of magnitude better times on "real" machines.