Finding connected components in 2D integer coordinate system

Question

Let L be a finite set of pair of integers.

I want to partition L according to connectivity

The connectivity is intuitive :

Rule1 (l1,l2),(l3,l4) is connected if (l1-l3, l2-l4) equals to one of (0,1),(0,-1),(1,0),(-1,0).

Rule2 If (l1,l2),(l3,l4) is connected and (l3,l4),(l5,l6) is connected, then (l1,l2),(l5,l6) is connected

(auxiliary rule : Only consider Rule1, Rule2.)

I know two methods, but I am not satisfied perfectly for both methods.

method 1)
MorphologicalComponents
As far as I know the method is a graphical approach. So to apply it to following example :

Connected components of {{1,1},{1,2},{1234512345,1234512345},{1234512345,1234512346}} ?

The answer is of course very easy :

Answer : {{{1,1},{1,2}},{{1234512345,1234512345},{1234512345,1234512346}}}

but to use MorphologicalComponents, we should introduce a rectangle of dimension 1234512345*1234512346. We start with so many unnecessary pixels. That's the point I am not satisfied.

method 2) ConnectedComponents

ConnectedComponents[{1 <-> 2, 1 <-> 3, 4 <-> 5}]
{{2, 1, 3}, {5, 4}}

To apply this method to solve my problem(mentioned on the beginning of this post), we should assert propositions such as followings, one by one. It will put a lot of strain on the RAM.

{50,50} and {50,51} is connected
{50,51} and {50,52} is connected
{50,52} and {50,53} is connected
...
{60,21} and {60,22} is connected
{60,22} and {60,23} is connected
...
{20,20} and {20,21} is connected
{20,21} and {20,22} is connected
...

Also the kernel of mathematica would look for this Connectivity Table everytime to get the job done. There is no way to explain mathematica kernel what is connectivity?, other then offering this bunch of table.

So I need yet another method. Can you suggest a good method ? I am almost completely new to this kind of field (connectivity, graph theory, MorphologicalComponents)

If possible, you can also answer for the case when diagonal connectivity is also allowed. Thank you!

+-+-+-+-+-+-+- After Henrik Schumacher's answer +-+-+-+-+-+-+-

It works well. Thank you!

You can get the picture of Euglena in mathematica help file ref/ MorphologicalComponents.

But one a little sad thing is it is fairly slower then MorphologicalComponents. For example,

MorphologicalComponents takes 0. seconds
Schumacher's method takes 0.15625 seconds

I will update the post soon!

Answer 1

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...