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 clustering[config_] := Module[{output, csizes, cindices}, output = MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False]; csizes = Rest@Sort@Tally@Flatten@output; cindices = {First@#, Position[output, First@#]} & /@ csizes; {csizes, cindices, output} ]; and apply: inputConfig = {{-1, -1, -1, 0}, {-1, 0, 0, -1}, {-1, -1, 0, -1}, {0, -1, 0, 0}}; clustering@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.