# Finding connected components in an array

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

• Oh, and it feels a bit overkill/inefficient(?) to use connectedcomponents for graphs, Commented Feb 15, 2014 at 16:09
• I'm not sure I understand the connectedness criterion completely. Are you looking for MorphologicalComponents? The same with a custom neighbourhood definition? You say "if they either are on top of another, or one is down-left of the other" but this is not what your code does. Please make a sketch or ensure that the description is very clear (i.e. "these are neighbours, those are not neighbours). Commented Feb 15, 2014 at 16:16
• Ok, added a better description i think. Commented Feb 15, 2014 at 16:40
• So yeah, it is sort of Morphological components, if you tilt your head 45 degrees left. Commented Feb 15, 2014 at 17:33
• In the second paragraph you wrote: "Two coordinates are neighbouring if they either are on top of another, or one is down-left of the other." Is the "on top of another" part a mistake? "On top" would suggest to me that coordinates in alternate rows (like the first and third row) which appear to be along the same vertical line in the parallelogram diagram are neighbors..
– Aky
Commented Feb 16, 2014 at 7:38

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


• +1 for idea - I did similar, like yours was tragically slow on other than toy-sized inputs, mulling how to do this efficiently.
– ciao
Commented Feb 16, 2014 at 4:41
• Looks nice! I can't really follow the code, but how do you make the output use the coordinates in l1? Your output also give coordinates to points which do not correspond to elements in l1... Is it possible to work in the native coordinates of l1, or should one simply "fix" your output somehow? Commented Feb 16, 2014 at 9:10
• I guess it suffices to add clusters = clusters /. {{r_Integer, c_Integer} :> If[OddQ[r + c], Sequence @@ {}, {r, (c - r)/2 + 1}]}; Commented Feb 16, 2014 at 9:29
• Ok, I modified it a bit, see the code in the bottom of the post. Thank you! Works like a charm! Commented Feb 16, 2014 at 10:02
• @Paxinum Glad you could use it. I made it with a few pieces I had already done for solving other problems:) Commented Feb 16, 2014 at 17:10

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").

• Looks interesting, I will try it out tomorrow :). Good performance there! Commented Feb 20, 2014 at 19:00
• @Paxinum: I hope you find it interesting, and I'm curious re: timings in your environment. Hold on to your hat - I just cobbled up a new method that looks to be over 3X faster, shorter and cleaner code, much lower resource needs, and I think I can double performance with some tweaks - very late here, so will be in spare time later today...
– ciao
Commented Feb 21, 2014 at 11:11
• Oh, I am looking forward to that post! Commented Feb 21, 2014 at 14:40
• @Paxinum: Updated. Pretty happy with results...
– ciao
Commented Feb 22, 2014 at 0:51
• @Paxinum: Well, that was short-lived. Replaced with stupid-fast update.
– ciao
Commented Feb 22, 2014 at 11:10

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.

• Notice that I modified the accepted answer a bit. Can you please also account for turning your output into the specified form? I tied your code, and it is slower on my machine, (but then, I added code for changing the output.) Commented Feb 16, 2014 at 11:41
• I added onepts = #[[3]] & /@ weirdNeighbors[array]; tiles = First[#[[4 ;;]]] & /@ weirdNeighbors[array]; tiles = Join @@ tiles; onepts = List /@ (Join @@ onepts); Join[onepts, tiles]; Commented Feb 16, 2014 at 11:41
• But, the arrays i am interested in, as in the example, all rows are weakly decreasing, and also the same holds for the two diagonals (in some directions)... this might affect which algorithm performs better... Commented Feb 16, 2014 at 11:45
• @Paxinum: Feel free to do with is as you please, I'm not planning on updating it (or working further on a sparse-array method 20X faster yet), anything posted is public-domain. No idea why your timings differ, my timings in other queries usually fairly consistent with "real" machines I use and those of others - I would be curious to see your results with the example RandomInteger[{0, 100}, {200, 200}] test in your environment.
– ciao
Commented Feb 16, 2014 at 11:46