10
$\begingroup$

Here is how to Generate data

the function is below with two images to generate data

img1

enter image description here

img2 enter image description here

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):

enter image description here

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

Improve this question
$\endgroup$
11
  • $\begingroup$ In the image attached the value 267 in the first line delineates that the object (label 1) in seg1 intersects an object labeled 5 in seg2 at 267 positions $\endgroup$
    – Ali Hashmi
    May 21, 2017 at 22:36
  • $\begingroup$ So are the two 2d arrays seg1 and seg2? In your pos printout, if each column reflects a single integer label and each row also reflects a unique integer label, how can you have two nonzero values in the first column, 204 and 118? $\endgroup$
    – MikeY
    May 21, 2017 at 22:37
  • $\begingroup$ @MikeY sorry i read your comment late. The thing is the labeled data in seg1 is 1,2,3....51 and seg2 is labeled with 1,2,3....50 so the entries being non-zero means that the 4th and 6th labels in seg1 intersect with the label 1 in seg2 $\endgroup$
    – Ali Hashmi
    May 21, 2017 at 22:56
  • $\begingroup$ maybe SparseArray[SparseArray@Normal[Rest@Counts[Flatten[(1 - Unitize[seg1 - #]) seg2]]]&/@Range[51]]// Normal? $\endgroup$
    – kglr
    May 21, 2017 at 23:24
  • $\begingroup$ @kglr let me add the raw images and the function from which the arrays are created. because currently your suggestion is giving a few errors. $\endgroup$
    – Ali Hashmi
    May 21, 2017 at 23:28

7 Answers 7

Reset to default
7
$\begingroup$

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

Improve this answer
$\endgroup$
8
  • $\begingroup$ @AliHashmi It works fine for me. Perhaps you need to restart the kernel or clear some variables? $\endgroup$
    – Carl Woll
    May 22, 2017 at 14:35
  • $\begingroup$ your method works like a charm !! plus one'd it already :) $\endgroup$
    – Ali Hashmi
    May 22, 2017 at 16:27
  • $\begingroup$ @AliHashmi IMHO you really should change the Accept to this post. I somehow overlooked it completely before starting on an answer myself, and it seems to be far and away the best method provided. $\endgroup$
    – Mr.Wizard
    May 22, 2017 at 16:32
  • $\begingroup$ @CarlWoll there is a slight problem with your approach. Your method here gives a dimension of {50,50} whereas it should have been {51,50} $\endgroup$
    – Ali Hashmi
    Jan 22, 2018 at 11:11
  • 1
    $\begingroup$ @AliHasmi The dimensions are {50, 50} because there is no intersection of the 51 from seg1 with any of the integers from seg2. This could be fixed by adding a dimension specification to the SparseArray code. That is, using SparseArray[Rule @@@ t, {Max[r1], Max[r2]}] instead of SparseArray[Rule @@@ t]. $\endgroup$
    – Carl Woll
    Jan 22, 2018 at 19:19
9
$\begingroup$ 
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

Improve this answer
$\endgroup$
3
  • $\begingroup$ this is brilliant stuff ! $\endgroup$
    – Ali Hashmi
    May 21, 2017 at 23:57
  • 2
    $\begingroup$ Thank you @Ali. You might want to wait a while before accepting as questions with accepted answers are less likely to attract attention. $\endgroup$
    – kglr
    May 22, 2017 at 0:28
  • $\begingroup$ i wish i could go back and double like this one :) $\endgroup$
    – Ali Hashmi
    May 22, 2017 at 11:30
5
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ nice work with PositionIndex. +1 $\endgroup$
    – Ali Hashmi
    May 22, 2017 at 11:08
5
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ This is still about 3-4 times slower than my method. Also, I would pull the definition of pos outside of fn2, as there is no reason to define pos every time you run fn2. $\endgroup$
    – Carl Woll
    May 22, 2017 at 16:21
  • $\begingroup$ @Mr.Wizard this is great ! +1 $\endgroup$
    – Ali Hashmi
    May 22, 2017 at 16:24
  • $\begingroup$ @CarlWoll That's weird, I completely overlooked your post. Sorry! I agree about externalizing pos but I wanted self-contained functions for the sake of this post. $\endgroup$
    – Mr.Wizard
    May 22, 2017 at 16:29
4
$\begingroup$

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 *)
Improve this answer
$\endgroup$
5
  • $\begingroup$ Nice method. I believe I can make your code cleaner while retaining the same performance and output by using KeySelect. May I edit your post? $\endgroup$
    – Mr.Wizard
    May 22, 2017 at 16:09
  • $\begingroup$ You bet, it was my first foray into Key[ ] type functions. I tried to use KeySelect[ ], but struggled with keys that are lists. $\endgroup$
    – MikeY
    May 22, 2017 at 16:11
  • $\begingroup$ Please check my edit. $\endgroup$
    – Mr.Wizard
    May 22, 2017 at 16:13
  • $\begingroup$ Nice simplification, thanks! $\endgroup$
    – MikeY
    May 22, 2017 at 16:15
  • $\begingroup$ @MikeY +1 for your fast approach $\endgroup$
    – Ali Hashmi
    May 22, 2017 at 16:31
1
$\begingroup$

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} *)

Improve this answer
$\endgroup$
0
$\begingroup$

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]} 
]]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.