Does this do what you have in mind? pts = {{1, 1}, {1, 2}, {1234512345, 1234512345}, {1234512345, 1234512346}}; edges = Join @@ MapIndexed[ Thread[{#1, #2[[1]]}] &, Nearest[ pts -> "Index", pts, {\[Infinity], 1}, DistanceFunction -> ManhattanDistance ][[All, 2 ;;]] ]; SparseArray`StronglyConnectedComponents[SparseArray[edges ->_]] For diagonal connectivity the same code with `ManhattanDistance` replaced by `ChessboardDistance` should work. **Edit** At least for the nondiagonal connectivity, this can be sped up by a factor of ten by avoiding `Nearest`. We can do so because we can find out the horizontal neighbors by going through the rows of sparse image matrix. This is easy to do for a matrix in CSR (compressed sparse row) format. For finding out the vertical neighbors, we just do that for the transposed image, too. First a couple of helper functions: QuickSparseArray[{rp_?VectorQ, ci_?MatrixQ, vals_?VectorQ, dims_?VectorQ, background_ : 0}] := With[{data = {Automatic, dims, background, {1, {rp, ci}, vals}}}, SparseArray @@ data ]; ThreadCount[] := "ParallelThreadNumber" /. ("ParallelOptions" /. SystemOptions["ParallelOptions"]); JobPointers[jobCount_Integer?Positive, threadCount_Integer?Positive] := Ceiling[Subdivide[0, jobCount, Min[threadCount, jobCount]]]; cFindHorizontalNeighbors = Compile[{{rp, _Integer, 1}, {ci, _Integer, 2}, {idx, _Integer, 1}, {start, _Integer}, {end, _Integer}}, Block[{bag, col, nextcol, i, j}, bag = Internal`Bag[Most[{0}]]; Do[ Do[ col = Compile`GetElement[ci, k, 1]; nextcol = Compile`GetElement[ci, k + 1, 1]; If[nextcol == col + 1, i = Compile`GetElement[idx, k]; j = Compile`GetElement[idx, k + 1]; Internal`StuffBag[bag, i]; Internal`StuffBag[bag, j]; Internal`StuffBag[bag, j]; Internal`StuffBag[bag, i]; ] , {k, Compile`GetElement[rp, row] + 1, Compile`GetElement[rp, row + 1] - 1}] , {row, start, end}]; Partition[Internal`BagPart[bag, All], 2] ], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True, RuntimeOptions -> "Speed" ]; Now the actual computations: n = Length[A["NonzeroValues"]]; B = QuickSparseArray[A["RowPointers"], A["ColumnIndices"], Range[n], Dimensions[A]]; BT = Transpose[B]; Bptr = JobPointers[Length[B], ThreadCount[]]; BTptr = JobPointers[Length[BT], ThreadCount[]]; edges = Join[ Join @@ cFindHorizontalNeighbors[B["RowPointers"], B["ColumnIndices"], B["NonzeroValues"], Most[Bptr] + 1, Rest[Bptr]], Join @@ cFindHorizontalNeighbors[BT["RowPointers"], BT["ColumnIndices"], BT["NonzeroValues"], Most[BTptr] + 1, Rest[BTptr]] ]; components = SparseArray`StronglyConnectedComponents[SparseArray[edges -> _, {n, n}]]; Bonus: We can get the colored sparse array with colors = Normal[SparseArray[ Join @@ components -> Join @@ (Range[Length[components]] Unitize[components]), n]]; Acolored = QuickSparseArray[A["RowPointers"], A["ColumnIndices"], colors, Dimensions[A]]; Now we can execute `Colorize[Acolored]` and behold. In principal, one could do a similar thing for the diagonal lookup. I am just not in the modd for that...