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:
Build an adjacency matrix for the colors that touch modulo periodicity.
Compute the connected components of this matrix with SparseArray`StronglyConnectedComponents (usually considerably faster than using Graph tools such as ConnectedComponents).
Create a lookup table cols in which cols[[c+1]] contains the new integer index for the old integer index c.
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
]
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
]
Notice the option setting CornerNeighbors -> False in periodicBoundaryCorrection; with the default CornerNeighbors -> True we would obtain
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]
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]
];
MorphologicalComponentsdoes take aPaddingoption, butPadding -> "Periodic"does not work. $\endgroup$