Skip to main content
3 of 3
deleted 2 characters in body
Henrik Schumacher
  • 109.5k
  • 7
  • 186
  • 323

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 principle, one could do a similar thing for the diagonal lookup. I am just not in the modd for that...

Henrik Schumacher
  • 109.5k
  • 7
  • 186
  • 323