intersection between two 2D arrays with labeled data is slow

Question

Here is how to Generate data

the function is below with two images to generate data

img1

img2

segmentImage[binarizedMask_?ImageQ] := 
Module[{seg, areas, indexMaxarea, maxArea},
seg = MorphologicalComponents@*ColorNegate@Dilation[binarizedMask, 1];
areas = ComponentMeasurements[seg, "Area"];
{indexMaxarea, maxArea} = First@MaximalBy[areas, Last] /. Rule -> List;
If[maxArea > 20000, ArrayComponents[seg, Length@areas, indexMaxarea -> 0], seg] ~
Dilation~1(* 20000 is just a threshold *)
];

generating data

{seg1, seg2} = segmentImage /@ {img1, img2};

Now the question

I have two 2-D arrays with dimensions {651,823}. Each array has labelled components i.e. at each position the value is either zero or an integer label. I have 51 objects in the first seg1 and 50 objects in the second array seg2.

I wish to calculate the number of positions where a given object from the first array intersects the objects in the second array.

I have concocted this strategy but it seems terribly slow:

(pos = Table[
 Length@Intersection[Position[seg1, i], Position[seg2, j]], {i, 
  51}, {j, 50}];) // AbsoluteTiming
(* {143.607, Null} *)

Since these arrays need to be modified and consequently the intersections need to recalculated in an iteration later, I find that this is not the best way to do it.

Any ideas of speeding the process up will be welcomed.

the output i have looks something like this (each row represents an integer label in seg1 and each column represents a label in seg2 and the values represent the number of positions the labelled objects in both arrays intersect):

Note: In the image attached the values 267 and 564 in the first line delineates that the object (label 1) in seg1 intersects an object labeled 5 and 7 in seg2 at 267 and 564 positions respectively

Answer 1

AbsoluteTiming[
 counts = Normal[SparseArray[SparseArray[#, {50}] &@ DeleteCases[
         Normal[Counts[Flatten[(1 - Unitize[seg1 - #]) seg2]]], 0 -> _] & /@ Range[51]]];]

{0.30279, Null}

(pos = Table[Length@Intersection[Position[seg1, i], Position[seg2, j]], {i, 
        51}, {j, 50}];) // AbsoluteTiming

{122.747, Null}

counts == pos

True

Answer 2

Something like the following would be pretty fast:

intersections[r1_, r2_] := Module[{t},
  t = DeleteCases[Tally@Flatten[Unitize[r1 r2] (100 r1 + r2)], {0, _}];
  t[[All, 1]] = IntegerDigits[t[[All, 1]], 100];
  Normal[SparseArray[Rule @@@ t]]
]

This assumes that the labels are all less than 100. For the OP example:

res=intersections[seg1, seg2]; //AbsoluteTiming
res[[;;5, ;;5]]

{0.006349, Null}

{{0, 0, 0, 0, 267}, {0, 0, 6888, 0, 0}, {0, 6511, 0, 0, 0}, {204, 0, 0, 7249, 70}, {0, 0, 0, 110, 3971}}

Answer 3

A simpler solution using PositionIndex

a1 = Flatten[seg1] // PositionIndex;
a2 = Flatten[seg2] // PositionIndex;
a3 = Table[Length@Intersection[a1[i], a2[j]], {i, 51}, {j, 50}]; // AbsoluteTiming

{0.930523, Null}

counts == pos == a3

True

Can't beat kglr's solution, though. It takes 0.41 s.

Answer 4

Pick is optimized to work on packed arrays. Let's try it. As a self-contained function:

fn[x_, y_] :=
 Module[{a, b, m, n, toRow},
   {m, n} = Max /@ {x, y};
   {a, b} = Flatten /@ {x, y};
   toRow = SparseArray[Normal @ KeyDrop[0] @ Counts @ #, {n}] &;
   Array[toRow @ Pick[b, a, #] &, m]
 ]

Output matches kglr's:

img1 = Import["https://i.stack.imgur.com/oKpri.png"];
img2 = Import["https://i.stack.imgur.com/yzvfv.png"];

{seg1, seg2} = segmentImage /@ {img1, img2};

fn[seg1, seg2] == counts  (* True *)

Performance is superior by more than an order of magnitude:

fn[seg1, seg2] // RepeatedTiming // First

Normal[SparseArray[
    SparseArray[#, {50}] &@
       DeleteCases[Normal[Counts[Flatten[(1 - Unitize[seg1 - #]) seg2]]], 
        0 -> _] & /@ Range[51]]] // RepeatedTiming // First

0.026
0.488

---

The code above still scans the data fifty times. That's usually a bad idea. With more segments this might become an issue. As an attempt to avoid the problem I'll use Szabolcs's function from How to efficiently find positions of duplicates? (as pos) as a substitute for Pick:

fn2[x_, y_] :=
  Module[{a, b, n, pos, toRow},
    n = Max[y];
    {a, b} = Flatten /@ {x, y};
    pos[lst_] := GatherBy[Range @ Length @ lst, lst[[#]] &];
    toRow = SparseArray[Normal @ KeyDrop[0] @ Counts @ #, {n}] &;
    toRow @ b[[#]] & /@ Rest @ pos @ a
  ]

Test:

fn[seg1, seg2] == fn2[seg1, seg2]

fn2[seg1, seg2] // RepeatedTiming // First

True
0.024

Only a slight edge over fn here but I expect it to pull ahead as the number of segments rises.

Answer 5

OK, as a learning exercise in associations I worked through this, and ended up with a faster result than the @kglr one. Count me surprised!

RepeatedTiming[
  cr = {seg1 // Flatten, seg2 // Flatten} // Transpose // Counts; 
  sa = SparseArray[ KeySelect[cr, FreeQ[0]] // Normal] // Normal;
]

(* out[11]={0.153, Null} *)

@kglr's

AbsoluteTiming[counts = Normal[SparseArray[SparseArray[#, {50}] &@DeleteCases[
   Normal[Counts[Flatten[(1 - Unitize[seg1 - #]) seg2]]], 0 -> _] & /@ Range[51]]];]

(* out[11]={0.478285, Null} *)

The matrices match (minus the fact that I did not pad the sa matrix... the last row is all zeros)

sa == Counts // Most
(* true *)

Answer 6

This approach works for the posted question as well as for cases when there are gaps in labels. What i mean by gaps is that the first or the second segmented mask may have a structure as follows: {1,2,3 ... 10,13,14,20} (11,12 as well as 16-19 are the missing labels).

overlapMatrix[seg1_,seg2_]:=Block[{keys1,keys2,mask,map,rules1,rules2},
keys1 = Keys@ComponentMeasurements[seg1, "Label"];
keys2 = Keys@ComponentMeasurements[seg2, "Label"];
mask= Unitize[seg1*seg2];
map = Normal@Counts@Thread[{SparseArray[mask*seg1]["NonzeroValues"],SparseArray[mask*seg2]["NonzeroValues"]}];
{rules1,rules2}= Dispatch@Thread[# -> Range@Length@#]&/@{keys1,keys2};
map[[All,1,1]]=map[[All,1,1]]/.rules1;
map[[All,1,2]]=map[[All,1,2]]/.rules2;
Normal@SparseArray[map,{Length@keys1,Length@keys2}]
];

Applying it over the above images

overlapMatrix[seg1,seg2]

(* {0.167852, Null} *)

Answer 7

Transpose accepts a second parameter, which you can use to shuffle the values around a bit. I have used it to merge the value of the two arrays together in a list.

Normal[SparseArray[
   DeleteCases[Normal[
       Counts[Flatten[Transpose[{seg1, seg2}, {3, 2, 1}], {1, 2}]]],
       Rule[{0,0}|{0,_}|{_,0}, _]
   ], 
   {Max[seg1], Max[seg2]} 
]]