Return to Revisions

Built-in option

This sidesteps most of your code, so it might not be what you are looking for, but I believe your goal can be achieved with Mathematica's built-in image processing capability, specifically: MorphologicalComponents!

Define a new clustering function

clustering1[config_] := Module[{output, csizes, cindices},
   output =  MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False];
   csizes = Rest@Sort@Tally@Flatten@output;
   cindices = Module[
      {sa = SparseArray[output], xx, yy, sa1, sa2},
      sa1 = sa["NonzeroValues"];
      xx = GatherBy[Range@Length@sa1, sa1[[#]] &];
      sa2 = sa["NonzeroPositions"];
      yy = sa2[[#]] & /@ xx;
      Transpose[{sa1[[xx[[All, 1]]]], yy}]];
     ];
   {csizes, cindices, output}
  ];

and apply:

inputConfig = {{-1, -1, -1, 0}, {-1, 0, 0, -1}, {-1, -1, 0, -1}, {0, -1, 0, 0}};
clustering1@inputConfig

{{{1, 7}, {2, 2}}, 

 {{1, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}}, {2, {{2, 4}, {3, 4}}}},
 {{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}

The first item in the output is the list of {clusterNumber, clusterSize}, the second the list of {clusterNumber, clusterIndices}, and the third is the input array with cluster numbers replacing filled sites.

EDIT

It appears that finding the cluster indices with Position is extremely slow with large datasets, although the rest of the code is relatively fast. With thanks to ciao (see below), I've replaced it with a much faster construction.

Mimicking the recursive algorithm

Here is my take on what the OP was actually asking for: a Mathematic approach to a recursive algorithm. The actual algorithm is not much changed - the use of Sow and Reap in the OP's code is a good idea, and I agree that the recursion needed is not amenable to Nest or Fold - all I can do is clean it up a little and exchange the Do for a Scan over the populated sites:

clustering2[config_] := Module[
   {output = config,
    cnum = 0, length = Length@config,
    csearch, clusters, c},
   
   csearch[{i_, j_}, cnum_] := If[
     output[[i, j]] == -1,
     output[[i, j]] = cnum;
     Sow[{i, j}];
     csearch[#, cnum] & /@ Select[
       {i, j} + # & /@ {{0, 1}, {0, -1}, {1, 0}, {-1, 0}},
       (1 <= First@# <= length && 1 <= Last@# <= length) &]
    ];
   
   clusters = Reap[
      Scan[
       If[output[[Sequence @@ #]] == -1,
          cnum++;
          c = Reap[csearch[#, cnum]][[2, 1]];
          Sow[{cnum, Length@c, c}]] &,
       SparseArray[config]["NonzeroPositions"]]
    ][[2, 1]];
   
   {clusters, output}
  ];

Applying this to inputConfig gives

{{{1, 7, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}}, 

  {2, 2, {{2, 4}, {3, 4}}}}, 
 {{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}

The first item in the output is the list of {clusterNumber, clusterSize, clusterIndices} and the second is the input array with cluster numbers replacing filled sites.

I've run this on a 1000x1000 array with filling factor ~0.5, and although it is about 3 times slower than the first method, it does work.

Comments

I'd be interested to know if anybody can think of a more clever way of implementing the algorithm. Mathematica has the SparseArray object, which is a good fit for the input and output configurations we are dealing with here, so one might consider using those from the beginning for the input. Thanks to ciao, who pointed out that one can extract the positions of non-zero elements easily from a SparseArray with ["NonzeroElements"], we can restrict our attention to those only.