14
$\begingroup$

I would like to find the connected(with periodic boundary conditions) components in a large binary matrix. What I've tried, and it does a fair job is:

MorphologicalComponents[RandomInteger[{0, 1}, {100, 100}],CornerNeighbors -> False]

For instance:

ArrayPlot[MorphologicalComponents[RandomInteger[{0, 1}, {20, 20}], CornerNeighbors -> False], ColorFunction -> "Rainbow", ImageSize -> Large]

enter image description here

The only option missing here is the periodic boundary condition. MorphologicalComponents doesn't have a DistanceFunction to impose periodic boundary conditions there.

The other idea I had was to give up on MorphologicalComponents and use ClusteringComponents with a complicated DistanceFunction, but that is just a pain in the neck.


Example:

Let me give an example of what I mean by periodic boundary condition. For instance given the input:

{{1, 0, 0, 1, 1}, {0, 0, 1, 0, 0}, {1, 1, 1, 0, 1}, {0, 1, 0, 0, 0}, {1, 0, 0, 1, 0}}

I would expect to get something like:

{{1, 0, 0, 1, 1}, {0, 0, 2, 0, 0}, {2, 2, 2, 0, 2}, {0, 2, 0, 0, 0}, {1, 0, 0, 1, 0}}

As an output.

The arrayplots for these two matrices are:

enter image description here

enter image description here

$\endgroup$
6
  • $\begingroup$ Can you show an example of what you would like the output to look like? (even if using a smaller grid). I'm having trouble visualizing what you have in mind. $\endgroup$
    – DavidC
    Commented Mar 29, 2014 at 2:20
  • $\begingroup$ @DavidCarraher Sure. Just give me few minutes. $\endgroup$
    – Ali
    Commented Mar 29, 2014 at 2:30
  • 2
    $\begingroup$ @David: I think he just wants the grid to be treated as periodic, so that a cell adjacent to the top edge and a cell adjacent to the bottom edge with the same horizontal coordinate should be connected (and similarly for the left and right edges). In this example, the five-cell dark blue component and the two-cell purple component in the top middle-left should both be parts of the large orange component at bottom left. $\endgroup$
    – user484
    Commented Mar 29, 2014 at 2:37
  • $\begingroup$ Ok. I now see it. I needed to remove the morphological components and look at the 1's and 0's to fully get the picture. $\endgroup$
    – DavidC
    Commented Mar 29, 2014 at 2:41
  • 2
    $\begingroup$ As it turns out, MorphologicalComponents does take a Padding option, but Padding -> "Periodic" does not work. $\endgroup$
    – user484
    Commented Mar 29, 2014 at 2:41

5 Answers 5

17
$\begingroup$

This generates a 20 by 20 binary matrix and finds the morphological components.

SeedRandom[11];
m=RandomInteger[{0,1},{20,20}];
a=MorphologicalComponents[m,CornerNeighbors->False]

grid

Notice that morphological component 2, in row 1, col 6, abuts morphological component 42 in row 20, col 6. Morphological components 2 and 39 abut in column 8.

These correspond to the abutments of dark blue (mc2) and red (mc42), at the top and bottom of column 6, and dark blue (mc2) and orange (mc39), in column 8, as shown in the ArrayPlot below without periodic boundary conditions imposed:

ArrayPlot[a,ColorFunction->"Rainbow",ImageSize->400]

mc1


The following will display the same binary matrix data with periodic boundary conditions.

periodicBoundaryArrayPlot[morphcomponents_]:=
Module[{connectedMorphologicalComponents},
    connectedMorphologicalComponents[mc_]:=Cases[Union@Partition[Riffle[mc[[1]],mc[[-1]]],2],({x_,y_}/;x!=0&&y!=0&&x!=y):> 
    UndirectedEdge@@Sort[{x,y}]];
    ArrayPlot[Replace[morphcomponents,Flatten[Thread[Rule[Most@#,Last@#]]&/@
    ConnectedComponents[Graph[Union[connectedMorphologicalComponents[morphcomponents],
    connectedMorphologicalComponents[Transpose[morphcomponents]]]]]],2],
    ColorFunction->"Rainbow",ImageSize->400]]

The following relies on same initial data as above, but now takes into account the periodic boundary conditions.

periodicBoundaryArrayPlot[a]

mc2


Analysis

In the original code the vertical and horizontal abutments were not explicitly named. They are named here to facilitate interpretation of the code.

connectedMorphologicalComponents[mc_] := 
 Cases[Union@
   Partition[Riffle[mc[[1]], mc[[-1]]], 
    2], ({x_, y_} /; x != 0 && y != 0 && x != y) :> 
       UndirectedEdge @@ Sort[{x, y}]]
verticalAbutments = connectedMorphologicalComponents[a]
horizontalAbutments = connectedMorphologicalComponents[Transpose[a]]

undirected edges

This shows that morphological component 1 runs into mc 41 vertically, and so. on. The symbol between 1 and 41 stands for an undirected edge of a graph.

Here the graph of the above edges is made and the connected components of that graph (not the morphological components) are isolated.

ConnectedComponents[Graph[Union[verticalAbutments,horizontalAbutments]]]

{{4, 3, 36, 40, 37}, {42, 2, 39}, {1, 41}, {14, 11}}

This means that morphological components 4, 3, 36, 40, and 37 are actually a single component (because these components run into each other either vertically or horizontally). They should thus be colored identically.

Components 42, 2, and 39 should be reduced to a single component with a single color.

Likewise for components 1 and 41; and for 14 and 11.

Replace[morphcomponents,Flatten[Thread[Rule[Most@#,Last@#]] carries out the required replacements and reductions of components. It replace the first components in a sublist with the last component in the sublist:

{4-> 37, 3-> 37, 36->37, 40->37, 42->39, 2->39, 1->41, 14->11}
$\endgroup$
3
  • 1
    $\begingroup$ This looks very nice! But it seems you're only considering periodicity in the vertical direction. Going by the second example in the question, the OP wants both the horizontal and the vertical boundaries of the image to be periodic. $\endgroup$
    – user484
    Commented Mar 29, 2014 at 9:37
  • $\begingroup$ The procedure now checks the horizontal and vertical borders. $\endgroup$
    – DavidC
    Commented Mar 29, 2014 at 13:28
  • $\begingroup$ @DavidCarraher Thanks for the very nice answer. $\endgroup$
    – Ali
    Commented Mar 29, 2014 at 21:51
7
$\begingroup$

I think bill's tiling and extracting idea in his deleted answer is actually nice, only a little more effort is needed.

First we define a handy plot function:

Clear[morphPlot]
morphPlot[m_] := ArrayPlot[m, ColorFunction -> "Pastel"] /. (List @@ ColorData["Pastel"][0]) -> {0, 0, 0}

We generate a test array m, tile it and apply the MorphologicalComponents:

m = RandomInteger[{0, 1}, {15, 20}];

largeMat = ArrayFlatten[{{0, m, 0}, {m, m, m}, {0, m, 0}}];
largeMorph = MorphologicalComponents[largeMat, CornerNeighbors -> False];

(* The replacement is for seperating similar colors from each other,
   not necessary for calculation: *)
largeMorph = largeMorph /. Dispatch[Thread[# -> RandomSample[#]] &@Range[Max@largeMorph]];

morphPlot[largeMorph]

tiled morphological components

then delete those components who don't have intersection with the central block (i.e. the original m):

morphPart = Partition[largeMorph, Dimensions[m]];
allLabels = Union[Flatten[largeMorph]];
trueLabels = Union[Flatten[morphPart[[2, 2]]]];
trueMorphPart = morphPart /. Dispatch[Thread[Complement[allLabels, trueLabels] -> 0]];

relevant components

The rest work is some replacements inside the equivalent classes:

(
 trueMorphPart = trueMorphPart /.
      Dispatch[
       MapThread[
                 If[#2 != 0, #1 -> #2, {}] &,
                 {trueMorphPart[[2, 2]], trueMorphPart[[##]] & @@ #},
                 2] // Flatten // Union
                ]
) & /@ {{1, 2}, {2, 1}, {2, 3}, {3, 2}};

periodicMorph = trueMorphPart[[2, 2]];

Show[{morphPlot[periodicMorph],
      MapIndexed[Text[Style[#1, 13], {#2[[2]], 15 - #2[[1]] + 1} - .5] &,
                 periodicMorph, {2}] // Flatten // Graphics
     }]

morphological components with periodic boundary condition

$\endgroup$
7
$\begingroup$

I wrote such a postprocessor for MorphologicalComponents as part of some other post and thought it might be a good idea to post it also here. The following looks sluggish but I put some effort into making it efficient. The algorithm consists of four steps:

  1. Build an adjacency matrix for the colors that touch modulo periodicity.

  2. Compute the connected components of this matrix with SparseArray`StronglyConnectedComponents (usually considerably faster than using Graph tools such as ConnectedComponents).

  3. Create a lookup table cols in which cols[[c+1]] contains the new integer index for the old integer index c.

  4. Use cols to lookup the new indices in parallelized way.

This is the code:

ClearAll[periodicBoundaryCorrection];
periodicBoundaryCorrection[
 A_?MatrixQ, 
 OptionsPattern[{CornerNeighbors -> True}]
 ] := Module[{maxA, a, b, c, d, pos, r1, r2, r3, r4, r5, r6, edges, α, β, colorcomp, cols},
   maxA = Max[A];
   
   a = A[[1]];
   b = A[[-1]];
   pos = DeleteCases[Range[Length[a]] Unitize[a b], 0];
   r1 = Sort /@ Transpose[{a[[pos]], b[[pos]]}];
   
   c = A[[All, 1]];
   d = A[[All, -1]];
   pos = DeleteCases[Range[Length[c]] Unitize[c d], 0];
   r2 = Sort /@ Transpose[{c[[pos]], d[[pos]]}];
   
   edges = Union[r1, r2];

   If[OptionValue[CornerNeighbors],

    pos = DeleteCases[Range[Length[a] - 1] Unitize[Rest[a] Most[b]], 0];
    r1 = Sort /@ Transpose[{a[[pos + 1]], b[[pos]]}];
    pos = DeleteCases[Range[Length[b] - 1] Unitize[Rest[b] Most[a]], 0];
    r2 = Sort /@ Transpose[{b[[pos + 1]], a[[pos]]}];

    pos = DeleteCases[Range[Length[c] - 1] Unitize[Rest[c] Most[d]], 0];
    r3 = Sort /@ Transpose[{c[[pos + 1]], d[[pos]]}];
    pos = DeleteCases[Range[Length[d] - 1] Unitize[Rest[d] Most[c]], 0];
    r4 = Sort /@ Transpose[{d[[pos + 1]], c[[pos]]}];
    
    α = c[[1]];
    β = d[[-1]];
    r5 = If[α β != 0, {Sort[{α, β}]}, {}];
    α = c[[-1]];
    β = d[[1]];
    r6 = If[α β != 0, {Sort[{α, β}]}, {}];
    
    edges = Union[edges, r1, r2, r3, r4, r5, r6]
   ];
   edges++;
   If[Length[edges] == 0,
    A,
    colorcomp = SparseArray`StronglyConnectedComponents[
      SparseArray[
       Join[edges, Transpose[Transpose[edges][[{2, 1}]]]] -> 
        1, {maxA + 1, maxA + 1}, 0]
      ];
    cols = Compile[{{idx, _Integer, 1}, {acc, _Integer, 1}},
       Block[{colors, j, threshold},
        colors = Table[0, {i, 1, acc[[-1]]}];
        j = 0;
        threshold = Compile`GetElement[acc, j + 1];
        colors[[idx]] = Table[
          If[i > threshold,
           j++;
           threshold = Compile`GetElement[acc, j + 1];
           ];
          j, {i, 1, Length[idx]}];
        colors
        ]
       ][
      Join @@ colorcomp,
      Accumulate[Length /@ colorcomp]
      ];
    Compile[{{a, _Integer, 1}, {cols, _Integer, 1}},
      cols[[a + 1]],
      RuntimeAttributes -> {Listable},
      Parallelization -> True
      ][A, cols]
    ]
   ];

Here's a usage example:

img = ExampleData[{"ColorTexture", "FloralPattern2"}];
GraphicsRow[
 {
  img,
  Colorize[
   MorphologicalComponents[img, 0.75]
   ],
  Colorize[
   periodicBoundaryCorrection[MorphologicalComponents[img, 0.75]]
   ]
  }, ImageSize -> 600
 ]

enter image description here

And here with random patterns:

SeedRandom[666];
A = RandomInteger[{0, 1}, {30, 30}];
GraphicsRow[{
  ArrayPlot[A],
  Colorize[
   MorphologicalComponents[A, CornerNeighbors -> False]
   ],
  Colorize[
   periodicBoundaryCorrection[
    MorphologicalComponents[A, CornerNeighbors -> False],
    CornerNeighbors -> False
    ]
   ]
  },
 ImageSize -> 600
 ]

enter image description here

Notice the option setting CornerNeighbors -> False in periodicBoundaryCorrection; with the default CornerNeighbors -> True we would obtain

enter image description here

Edit: MorphologicalEulerNumber

Not sure whether this works correctly. Consider it as a first shot:

periodicMorphologicalEulerNumber[A_?MatrixQ, 
   OptionsPattern[{CornerNeighbors -> True}]] := Subtract[
   Max[
    periodicBoundaryCorrection[
     MorphologicalComponents[A, 
      CornerNeighbors -> OptionValue[CornerNeighbors]], 
     CornerNeighbors -> OptionValue[CornerNeighbors]
     ]],
   Ramp[Max[
     periodicBoundaryCorrection[
      MorphologicalComponents[1 - A, 
       CornerNeighbors -> OptionValue[CornerNeighbors]], 
      CornerNeighbors -> OptionValue[CornerNeighbors]
      ]] - 1]
   ];

Some simple test:

n = 20;
A = DiskMatrix[2 n, 8 n];
A[[2 n + 1 ;; 6 n, 2 n + 1 ;; 6 n]] -= DiskMatrix[n, 4 n];
A = Join[
   A[[Quotient[Length[A], 2] + 1 ;;]],
   A[[;; Quotient[Length[A], 2]]]
   ];

Image[A]

MorphologicalEulerNumber[Image[A]]
periodicMorphologicalEulerNumber[A]

enter image description here

2

0

Edit: Attempt for 3D case

Follows the logic of the above, but it is totally untested:

periodicBoundaryCorrection[A_?(TensorQ[#] && TensorRank[#] == 3 &),
  OptionsPattern[{CornerNeighbors -> True}]] := 
 Module[{maxA, a, b, c, d, e, f, pos, r1, r2, r3, r4, r5, r6, 
   edges, \[Alpha], \[Beta], colorcomp, cols},
  maxA = Max[A];
  a = A[[1]];
  b = A[[-1]];
  pos = DeleteCases[Range[Length[a]] Unitize[a b], 0];
  r1 = Sort /@ Transpose[{a[[pos]], b[[pos]]}];
  c = A[[All, 1]];
  d = A[[All, -1]];
  pos = DeleteCases[Range[Length[c]] Unitize[c d], 0];
  r2 = Sort /@ Transpose[{c[[pos]], d[[pos]]}];
  e = A[[All, All, 1]];
  f = A[[All, All, -1]];
  pos = DeleteCases[Range[Length[e]] Unitize[e f], 0];
  r3 = Sort /@ Transpose[{e[[pos]], f[[pos]]}];
  edges = Union[r1, r2, r3];
  If[OptionValue[CornerNeighbors],
   Print["The option CornerNeighbors\[Rule]True is currently not implemented. Using CornerNeighbors\[Rule]False"];
   ];
  
  edges++;
  If[Length[edges] == 0, A, 
   colorcomp = 
    SparseArray`StronglyConnectedComponents[
     SparseArray[
      Join[edges, Transpose[Transpose[edges][[{2, 1}]]]] -> 
       1, {maxA + 1, maxA + 1}, 0]];
   cols = 
    Compile[{{idx, _Integer, 1}, {acc, _Integer, 1}}, 
      Block[{colors, j, threshold}, 
       colors = Table[0, {i, 1, acc[[-1]]}];
       j = 0;
       threshold = Compile`GetElement[acc, j + 1];
       colors[[idx]] = Table[If[i > threshold, j++;
          threshold = Compile`GetElement[acc, j + 1];];
         j, {i, 1, Length[idx]}];
       colors]][Join @@ colorcomp, Accumulate[Length /@ colorcomp]];
   Compile[{{a, _Integer, 1}, {cols, _Integer, 1}}, cols[[a + 1]], 
     RuntimeAttributes -> {Listable}, Parallelization -> True][A,cols]]
  ]

periodicMorphologicalEulerNumber[A_?(TensorQ[#] && TensorRank[#] == 3 &), 
   OptionsPattern[{CornerNeighbors -> True}]] := Subtract[
   Max[
    periodicBoundaryCorrection[
     MorphologicalComponents[A, 
      CornerNeighbors -> OptionValue[CornerNeighbors]], 
     CornerNeighbors -> OptionValue[CornerNeighbors]
     ]],
   Ramp[Max[
     periodicBoundaryCorrection[
      MorphologicalComponents[1 - A, 
       CornerNeighbors -> OptionValue[CornerNeighbors]], 
      CornerNeighbors -> OptionValue[CornerNeighbors]
      ]] - 1]
   ];
$\endgroup$
8
  • $\begingroup$ Hi chris! Thank you for the flowers. This is quite some time ago and I am quite busy today and the next couple of days. I think the 3D case just requires to adjust the computation of the edges. For this one just has to scan the borders of the 3D tensor A, i.e., compare A[[1]] to A[[-1]]; A[[All,1]] to A[[All,-1]]; and ; A[[All,All,1]] to A[[All,All,-1]]. The case CornerNeighbors is certainly much more complicated than in 2D... $\endgroup$ Commented Oct 15, 2021 at 7:38
  • $\begingroup$ I am not sure what you mean by MorphologicalEulerComponent, though. $\endgroup$ Commented Oct 15, 2021 at 7:39
  • $\begingroup$ Great post. Would you happen to have a simplified version which simply computes the MorphologicalEulerNumber (sorry for the wrong wording) for periodic images? Do you have a 3D version as well? – $\endgroup$
    – chris
    Commented Oct 15, 2021 at 7:41
  • 1
    $\begingroup$ See my recent edits. That's best I can do for now. (Have to run...) $\endgroup$ Commented Oct 15, 2021 at 8:15
  • 1
    $\begingroup$ @chris Obviously, it isn't! ^^ Well, if you give me a working definition, then I will see what I can do about it. $\endgroup$ Commented Oct 15, 2021 at 10:18
4
$\begingroup$
SeedRandom[43];
m = RandomInteger[{0, 1}, {14, 12}];
m1 = ArrayPad[m, 1, "Periodic"];
db = Dimensions@m1;
m1[[1, 1]] = m1[[1, db[[2]]]] = m1[[db[[1]], db[[2]]]] = m1[[db[[1]], 1]] = 0
b = MorphologicalComponents[m1, CornerNeighbors -> False];
t[{x_, y_}] := Flatten[{{{#, 1}, {#, y}} & /@ Range@x, {{1, #}, {x, #}} & /@  Range@y}, 1]
k = b //.((Min@#:>Max@#) & /@({b[[Sequence @@ #[[1]]]],b[[Sequence @@ #[[2]]]]} & /@ t[db]));

Show result:

GraphicsRow[{ArrayPlot[ MorphologicalComponents[m, CornerNeighbors -> False], ColorFunction -> "BrightBands"], 
             ArrayPlot[k[[2 ;; -2, 2 ;; -2]], ColorFunction -> "BrightBands"]}]

Mathematica graphics

$\endgroup$
1
  • 2
    $\begingroup$ I found a possible failure. $\endgroup$
    – Silvia
    Commented Mar 29, 2014 at 16:39
0
$\begingroup$

I came across a similar problem for 3D images and found Henrik's original untested 3D code not working properly. With some tweak of Henrik's code, I made the 3D version:

periodicBoundaryCorrection3D[A_?(TensorQ[#] && TensorRank[#] == 3 &), 
  OptionsPattern[{CornerNeighbors -> True}]] := 
 Module[{maxA, a, b, c, d, e, f, pos, r1, r2, r3, r4, r5, r6, 
   edges, \[Alpha], \[Beta], colorcomp, cols}, maxA = Max[A];
  a = A[[1]];
  b = A[[-1]];
  pos = Flatten[
    DeleteCases[
     CoordinateBoundsArray[{{1, Length[a]}, {1, Length[b]}}]  Unitize[
       a  b], {0, 0}, Infinity], 1];
  r1 = Sort /@ 
    Transpose[{Part[a, Sequence @@ #] & /@ pos, 
      Part[b, Sequence @@ #] & /@ pos}];
  
  c = A[[All, 1]];
  d = A[[All, -1]];
  pos = Flatten[
    DeleteCases[
     CoordinateBoundsArray[{{1, Length[c]}, {1, Length[d]}}]  Unitize[
       c  d], {0, 0}, Infinity], 1];
  r2 = Sort /@ 
    Transpose[{Part[c, Sequence @@ #] & /@ pos, 
      Part[d, Sequence @@ #] & /@ pos}];
  
  e = A[[All, All, 1]];
  f = A[[All, All, -1]];
  pos = Flatten[
    DeleteCases[
     CoordinateBoundsArray[{{1, Length[e]}, {1, Length[f]}}]  Unitize[
       e  f], {0, 0}, Infinity], 1];
  r3 = Sort /@ 
    Transpose[{Part[e, Sequence @@ #] & /@ pos, 
      Part[f, Sequence @@ #] & /@ pos}];
  
  edges = Union[r1, r2, r3];
  
  If[OptionValue[CornerNeighbors], 
   Print["The option CornerNeighbors\[Rule]True is currently not \
implemented. Using CornerNeighbors\[Rule]False"];];
  edges++;
  If[Length[edges] == 0, A, 
   colorcomp = 
    SparseArray`StronglyConnectedComponents[
     SparseArray[
      Join[edges, Transpose[Transpose[edges][[{2, 1}]]]] -> 
       1, {maxA + 1, maxA + 1}, 0]];
   cols = 
    Compile[{{idx, _Integer, 1}, {acc, _Integer, 1}}, 
      Block[{colors, j, threshold}, 
       colors = Table[0, {i, 1, acc[[-1]]}];
       j = 0;
       threshold = Compile`GetElement[acc, j + 1];
       colors[[idx]] = Table[If[i > threshold, j++;
          threshold = Compile`GetElement[acc, j + 1];];
         j, {i, 1, Length[idx]}];
       colors]][Join @@ colorcomp, Accumulate[Length /@ colorcomp]];
   Compile[{{a, _Integer, 1}, {cols, _Integer, 1}}, cols[[a + 1]], 
     RuntimeAttributes -> {Listable}, Parallelization -> True][A, 
    cols]]]

One example of 3D polycrystalline structures (left: MorphologicalComponents; right: periodicBoundaryCorrection3D)

left: MorphologicalComponents; right: periodicBoundaryCorrection3D

$\endgroup$

Your Answer

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

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