**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][1] (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][1], 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. [1]: http://mathematica.stackexchange.com/users/11467/ciao