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 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: enter image description here

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. GT-pattern

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;
   ];
share|improve this question
    
Oh, and it feels a bit overkill/inefficient(?) to use connectedcomponents for graphs, –  Paxinum Feb 15 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). –  Szabolcs Feb 15 at 16:16
    
Ok, added a better description i think. –  Paxinum Feb 15 at 16:40
    
So yeah, it is sort of Morphological components, if you tilt your head 45 degrees left. –  Paxinum Feb 15 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 Feb 16 at 7:38
show 1 more comment

3 Answers 3

up vote 4 down vote accepted
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}]}]

Mathematica graphics

share|improve this answer
    
+1 for idea - I did similar, like yours was tragically slow on other than toy-sized inputs, mulling how to do this efficiently. –  rasher Feb 16 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? –  Paxinum Feb 16 at 9:10
1  
I guess it suffices to add clusters = clusters /. {{r_Integer, c_Integer} :> If[OddQ[r + c], Sequence @@ {}, {r, (c - r)/2 + 1}]}; –  Paxinum Feb 16 at 9:29
1  
Ok, I modified it a bit, see the code in the bottom of the post. Thank you! Works like a charm! –  Paxinum Feb 16 at 10:02
    
@Paxinum Glad you could use it. I made it with a few pieces I had already done for solving other problems:) –  belisarius Feb 16 at 17:10
add comment

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.

enter image description here

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:

enter image description here

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

share|improve this answer
    
Looks interesting, I will try it out tomorrow :). Good performance there! –  Paxinum Feb 20 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... –  rasher Feb 21 at 11:11
    
Oh, I am looking forward to that post! –  Paxinum Feb 21 at 14:40
    
@Paxinum: Updated. Pretty happy with results... –  rasher Feb 22 at 0:51
    
@Paxinum: Well, that was short-lived. Replaced with stupid-fast update. –  rasher Feb 22 at 11:10
add comment

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]

enter image description here

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.

share|improve this answer
    
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.) –  Paxinum Feb 16 at 11:41
    
I added onepts = #[[3]] & /@ weirdNeighbors[array]; tiles = First[#[[4 ;;]]] & /@ weirdNeighbors[array]; tiles = Join @@ tiles; onepts = List /@ (Join @@ onepts); Join[onepts, tiles]; –  Paxinum Feb 16 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... –  Paxinum Feb 16 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. –  rasher Feb 16 at 11:46
add comment

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.