3
$\begingroup$

I have a $w \times h$ array with either $-1$ or $+1$.

I wish to monitor the size of the largest connected clusters with specific patterns that occur during a Monte Carlo simulation, i.e.,

  • a collection of $+1$'s,
  • a checkerboard pattern of $-1$ and $+1$'s, and,
  • horizontal and vertical stripes of $-1$ and $+1$'s.

Moreover, periodic boundary conditions apply such that opposite boundaries are connected.


Example arrays (checkerboard, stripes) can be imported by,

Uncompress[Import["https://pastebin.com/raw/0GaUtiFy"]]
Uncompress[Import["https://pastebin.com/raw/uF5UagqU"]] 

Ideally I am looking for a fast solution, such that monitoring every iteration of the Monte Carlo simulation will not slow it down to a crawl.

The only idea I have had so far was to use ComponentMeasurements but was not able to specify the patterns. Any help would be greatly appreciated.

$\endgroup$
2
$\begingroup$

This might not be the fasted method but it will get you started:

Preparation: The working horse is getComponent while periodicBoundaryCorrection only accounts for the fact that MorphologicalComponents does not allow for periodic boundary conditions. See the also here for a more detailed description.

checkermask = Array[Mod[Plus[##], 2] &, {3, 3}];
hormask = Array[Mod[#1, 2] &, {3, 3}];
vermask = Array[Mod[#2, 2] &, {3, 3}];

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]
    ]
   ];

getComponent[A_?MatrixQ, mask_] := periodicBoundaryCorrection[
  MorphologicalComponents[
   With[{img = Binarize[Image[A], 0]},
    With[{
      a = ImageCorrelate[img, mask/Total[mask, 2], Padding -> "Periodic"],
      b = ImageCorrelate[img, (1 - mask)/Total[1 - mask, 2], 
        Padding -> "Periodic"]
      },
     With[{ab = ImageMultiply[a, b]},
      Binarize[ImageAdd[ImageSubtract[a, ab], ImageSubtract[b, ab]], 
       0.99]
      ]
     ]
    ]
   ]
  ]

Application:

A = Uncompress[Import["https://pastebin.com/raw/0GaUtiFy"]];
B = Uncompress[Import["https://pastebin.com/raw/uF5UagqU"]];

GraphicsGrid[{
  {
   Image[A],
   Colorize[getComponent[A, checkermask]],
   Colorize[getComponent[A, hormask]],
   Colorize[getComponent[A, vermask]]
   },
  {
   Image[B],
   Colorize[getComponent[B, checkermask]],
   Colorize[getComponent[B, hormask]],
   Colorize[getComponent[B, vermask]]
   }
  }, ImageSize -> Large]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you, this definitely got me started. Reverse engineering your solution I found that getComponent[A_, mask_] := With[{B = -A ListCorrelate[mask, A, {2, 2}]}, ImageData[Binarize[Image[B], 3]]] almost does the trick. The only thing left is to implement the boundary conditions $\endgroup$ – user19218 Apr 11 '18 at 19:30
  • $\begingroup$ Ah sorry, I missed that. You have to set Padding - > "Periodic" in ListCorrelate. At least that should make the input for MorphologicalComponents correct. MorphologicalComponents won't match opposing components correctly. It has a Padding options but that does not allow for "Periodic"... $\endgroup$ – Henrik Schumacher Apr 11 '18 at 19:32
  • $\begingroup$ No sorry I will clarify, the components should found with periodic boundary conditions, i.e. in the checkerboard example there should only be one cluster found, instead of four $\endgroup$ – user19218 Apr 11 '18 at 19:40
  • $\begingroup$ In any case this question already has a answer $\endgroup$ – user19218 Apr 11 '18 at 19:52
  • $\begingroup$ Thanks for the hint. It was too late; already added a fix. $\endgroup$ – Henrik Schumacher Apr 11 '18 at 20:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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