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.
I'm not sure how well this will handle vary large datasets.