# Morphological Components with periodic boundary conditions

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] 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:  • 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. – DavidC Mar 29 '14 at 2:20
• @DavidCarraher Sure. Just give me few minutes. – Ali Mar 29 '14 at 2:30
• @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. – Rahul Mar 29 '14 at 2:37
• 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. – DavidC Mar 29 '14 at 2:41
• As it turns out, MorphologicalComponents does take a Padding option, but Padding -> "Periodic" does not work. – Rahul Mar 29 '14 at 2:41

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

SeedRandom;
m=RandomInteger[{0,1},{20,20}];
a=MorphologicalComponents[m,CornerNeighbors->False] 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] The following will display the same binary matrix data with periodic boundary conditions.

periodicBoundaryArrayPlot[morphcomponents_]:=
Module[{connectedMorphologicalComponents},
connectedMorphologicalComponents[mc_]:=Cases[Union@Partition[Riffle[mc[],mc[[-1]]],2],({x_,y_}/;x!=0&&y!=0&&x!=y):>
UndirectedEdge@@Sort[{x,y}]];
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] ## 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[], mc[[-1]]],
2], ({x_, y_} /; x != 0 && y != 0 && x != y) :>
UndirectedEdge @@ Sort[{x, y}]]
verticalAbutments = connectedMorphologicalComponents[a]
horizontalAbutments = connectedMorphologicalComponents[Transpose[a]] 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}

• 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. – Rahul Mar 29 '14 at 9:37
• The procedure now checks the horizontal and vertical borders. – DavidC Mar 29 '14 at 13:28
• @DavidCarraher Thanks for the very nice answer. – Ali Mar 29 '14 at 21:51

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}


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] 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]]; The rest work is some replacements inside the equivalent classes:

(
trueMorphPart = trueMorphPart /.
Dispatch[
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[], 15 - #2[] + 1} - .5] &,
periodicMorph, {2}] // Flatten // Graphics
}] SeedRandom;
m = RandomInteger[{0, 1}, {14, 12}];
db = Dimensions@m1;
m1[[1, 1]] = m1[[1, db[]]] = m1[[db[], db[]]] = m1[[db[], 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 @@ #[]]],b[[Sequence @@ #[]]]} & /@ t[db]));


Show result:

GraphicsRow[{ArrayPlot[ MorphologicalComponents[m, CornerNeighbors -> False], ColorFunction -> "BrightBands"],
ArrayPlot[k[[2 ;; -2, 2 ;; -2]], ColorFunction -> "BrightBands"]}] • I found a possible failure. – Silvia Mar 29 '14 at 16:39

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 SparseArrayStronglyConnectedComponents (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[];
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[];
β = d[[-1]];
r5 = If[α β != 0, {Sort[{α, β}]}, {}];
α = c[[-1]];
β = d[];
r6 = If[α β != 0, {Sort[{α, β}]}, {}];

edges = Union[edges, r1, r2, r3, r4, r5, r6]
];
edges++;
If[Length[edges] == 0,
A,
colorcomp = SparseArrayStronglyConnectedComponents[
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 = CompileGetElement[acc, j + 1];
colors[[idx]] = Table[
If[i > threshold,
j++;
threshold = CompileGetElement[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
] And here with random patterns:

SeedRandom;
A = RandomInteger[{0, 1}, {30, 30}];
GraphicsRow[{
ArrayPlot[A],
Colorize[
MorphologicalComponents[A, CornerNeighbors -> False]
],
Colorize[
periodicBoundaryCorrection[
MorphologicalComponents[A, CornerNeighbors -> False],
CornerNeighbors -> False
]
]
},
ImageSize -> 600
] Notice the option setting CornerNeighbors -> False in periodicBoundaryCorrection; with the default CornerNeighbors -> True we would obtain 