How to get the class of neighbors in a 2D periodic matrix?

Question

I'm trying to analyze the spatial structure of some rectangular arrays by counting the number and class/type of neighboring points around every point in the arrays, with periodic boundaries.

For instance, suppose I have a $m$ x $n$ array with 5 different "types" at equal frequencies:

(*Spatial domain and random uniform frequency of the five "types"*)
nTypes = 5;
m = 50;
n = 100;
freq = ConstantArray[1/nTypes, nTypes];

(*Matrix with all the individuals,representing the types as distinct \
integers*)
myMatrix = RandomInteger[RandomChoice[freq -> Range[nTypes]], {m, n}];

We can visualize the matrix like this, just to see the distribution of colors (I add Black for the "0" values generated by RandomInteger):

colors = {0 -> Black, 1 -> Blue, 2 -> Red, 3 -> Green, 4 -> Orange, 
   5 -> Purple};
MatrixPlot[myMatrix, ColorRules -> colors]

The question is, how can I efficiently get a list of lists, where each sublist contains the "type" and neighborhood its corresponding {m,n} point in the matrix?

Ideally, this would be in a function like this:

myNeighborhoodsList=[myMatrix_,radius_,neighborhoodType_]

With output like this (or similar):

{{{neighbor1->5},...,{lastneighbor->2}},...,{{neighbor1->3},...,{lastneighbor->4}}}

Where each sublist is an association between the position of a neighbor and its "type".

Ideally, the function should take some arbitrary radius, and be able to use either Moore or von Neumann neighborhoods (that's why I included these as arguments).

Currently I'm trying to map every {m,n} element with its corresponding neighborhood. For instance, for the element {m,n} and a Moore neighborhood of radius 1, with {m=5, n=6}, I'd do:

element = {5, 6};
{1 -> myMatrix[[element[[1]] - 1, element[[2]] - 1]], 
 2 -> myMatrix[[element[[1]], element[[2]] - 1]],
 3 -> myMatrix[[element[[1]] + 1, element[[2]] - 1]],
 4 -> myMatrix[[element[[1]] - 1, element[[2]]]],
 5 -> myMatrix[[element[[1]] + 1, element[[2]]]], 
 6 -> myMatrix[[element[[1]] - 1, element[[2]] + 1]],
 7 -> myMatrix[[element[[1]], element[[2]] + 1]],
 8 -> myMatrix[[element[[1]] + 1, element[[2]] + 1]]}

Giving a list in the form:

{1 -> 3, 2 -> 3, 3 -> 0, 4 -> 2, 5 -> 2, 6 -> 0, 7 -> 2, 8 -> 0}

I'm having trouble trying to generalize this to an arbitrary radius, and a different neighborhood (von Neumann, etc). I'm wondering if it's possible to use some sort of convolution kernel applied to the entire matrix with periodic boundaries, which would allow us to define arbitrary neighborhoods (not just Moore/von Neumann).

Answer 1

How about this?

extractNeighborsNeumann[m_, r_] := {
  center -> Part[m, r + 1, r + 1],
  neighbors -> extractNeighbors[
    m,
    DiamondMatrix[r] - CenterArray[1, {2 r + 1, 2 r + 1}]
    ]
  }
extractNeighborsMoore[m_, r_] := {
  center -> Part[m, r + 1, r + 1],
  neighbors -> extractNeighbors[
    m,
    BoxMatrix[r] - CenterArray[1, {2 r + 1, 2 r + 1}]
    ]
  }
extractNeighbors[m_, mask_] := Pick[Flatten[m], Flatten[mask], 1]

partitionMap[f_, m_, r_] := Developer`PartitionMap[f, 
  m, {2 r + 1, 2 r + 1}, {1, 1}, {r + 1, r + 1}]
myNeighborhoodsList[m_, r_, "Neumann"] := 
 partitionMap[extractNeighborsNeumann[#, r] &, m, r]
myNeighborhoodsList[m_, r_, "Moore"] := 
 partitionMap[extractNeighborsMoore[#, r] &, m, r]

m = Partition[Range[9], 3];
myNeighborhoodsList[m, 1, "Moore"]

(* Out: {{{center -> 1, neighbors -> {7, 3, 2, 4}}, {center -> 2, 
   neighbors -> {8, 1, 3, 5}}, {center -> 3, 
   neighbors -> {9, 2, 1, 6}}}, {{center -> 4, 
   neighbors -> {1, 6, 5, 7}}, {center -> 5, 
   neighbors -> {2, 4, 6, 8}}, {center -> 6, 
   neighbors -> {3, 5, 4, 9}}}, {{center -> 7, 
   neighbors -> {4, 9, 8, 1}}, {center -> 8, 
   neighbors -> {5, 7, 9, 2}}, {center -> 9, 
   neighbors -> {6, 8, 7, 3}}}} *)

myNeighborhoodsList[m, 1, "Neumann"]

(* Out: {{{center -> 1, neighbors -> {9, 7, 8, 3, 2, 6, 4, 5}}, {center -> 2, 
   neighbors -> {7, 8, 9, 1, 3, 4, 5, 6}}, {center -> 3, 
   neighbors -> {8, 9, 7, 2, 1, 5, 6, 4}}}, {{center -> 4, 
   neighbors -> {3, 1, 2, 6, 5, 9, 7, 8}}, {center -> 5, 
   neighbors -> {1, 2, 3, 4, 6, 7, 8, 9}}, {center -> 6, 
   neighbors -> {2, 3, 1, 5, 4, 8, 9, 7}}}, {{center -> 7, 
   neighbors -> {6, 4, 5, 9, 8, 3, 1, 2}}, {center -> 8, 
   neighbors -> {4, 5, 6, 7, 9, 1, 2, 3}}, {center -> 9, 
   neighbors -> {5, 6, 4, 8, 7, 2, 3, 1}}}} *)

Developer`PartitionMap is not well documented, also the documentation says that it has been superseded by BlockMap. However, a more thorough documentation of its arguments is available under Partition. Why not use BlockMap? BlockMap unfortunately does not support some of the arguments that we need to rely on in order to take the periodic boundaries into account.

Neighborhoods

The neighborhoods are implemented by creating "masks" of the elements belonging to the neighborhoods. A mask is a matrix of ones and zeros, where ones are neighbors.

Von Neumann

The Von Neumann neighborhood is implemented using DiamondMatrix, and CenterArray is used to remove the center:

r = 3;
m = DiamondMatrix[r] - CenterArray[1, {2 r + 1, 2 r + 1}];
m // MatrixPlot

Moore

The Moore neighborhood is implemented using BoxMatrix, and CenterArray is used to remove the center:

r = 3;
m = BoxMatrix[r] - CenterArray[1, {2 r + 1, 2 r + 1}];
MatrixPlot[m]

Other neighborhoods

Other neighborhoods can be constructed by combining functions such as DiamondMatrix, BoxMatrix, CrossMatrix, DiskMatrix, CenterArray etc.

Extension to get the indices of the matrix elements

It is easy to make a version of Developer`PartitionMap that also gives the index of the element it is currently mapping over. I'll illustrate it for the Moore neighborhood. Here is the partitionMapIndexed function:

partitionMapIndexed[f_, m_, r_] := MapIndexed[f,
  Partition[m, {2 r + 1, 2 r + 1}, {1, 1}, {r + 1, r + 1}],
  {2}
  ]

Now, we might do the following:

SetAttributes[myMod, Listable]
myMod[v_, max_] := If[
  v >= 1, Mod[v, max, 1],
  max - Mod[Abs[v], max]
  ]

neighborIndices[centerIndex_, width_, r_, mask_] := Pick[Flatten[Table[
    myMod[centerIndex + {i, j}, width],
    {i, -r, r},
    {j, -r, r}
    ], 1], Flatten[mask], 1]

extractNeighborsMoore[m_, idx_, width_, r_] := {
  center -> Part[m, r + 1, r + 1],
  centerIndex -> idx,
  neighbors -> extractNeighbors[
    m,
    BoxMatrix[r] - CenterArray[1, {2 r + 1, 2 r + 1}]
    ],
  neighborsIndices -> neighborIndices[
    idx,
    width,
    r,
    BoxMatrix[r] - CenterArray[1, {2 r + 1, 2 r + 1}]
    ]
  }

myNeighborhoodsList[m_, r_, "Moore"] := 
 partitionMapIndexed[extractNeighborsMoore[#, #2, Length[m], r] &, m, r]

m = Partition[Range[9], 3];
myNeighborhoodsList[m, 1, "Moore"]

(* Out: {{{center -> 1, centerIndex -> {1, 1}, 
   neighbors -> {9, 7, 8, 3, 2, 6, 4, 5}, 
   neighborsIndices -> {{3, 3}, {3, 1}, {3, 2}, {1, 3}, {1, 2}, {2, 
      3}, {2, 1}, {2, 2}}}, {center -> 2, centerIndex -> {1, 2}, 
   neighbors -> {7, 8, 9, 1, 3, 4, 5, 6}, 
   neighborsIndices -> {{3, 1}, {3, 2}, {3, 3}, {1, 1}, {1, 3}, {2, 
      1}, {2, 2}, {2, 3}}}, {center -> 3, centerIndex -> {1, 3}, 
   neighbors -> {8, 9, 7, 2, 1, 5, 6, 4}, 
   neighborsIndices -> {{3, 2}, {3, 3}, {3, 1}, {1, 2}, {1, 1}, {2, 
      2}, {2, 3}, {2, 1}}}}, {{center -> 4, centerIndex -> {2, 1}, 
   neighbors -> {3, 1, 2, 6, 5, 9, 7, 8}, 
   neighborsIndices -> {{1, 3}, {1, 1}, {1, 2}, {2, 3}, {2, 2}, {3, 
      3}, {3, 1}, {3, 2}}}, {center -> 5, centerIndex -> {2, 2}, 
   neighbors -> {1, 2, 3, 4, 6, 7, 8, 9}, 
   neighborsIndices -> {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 3}, {3, 
      1}, {3, 2}, {3, 3}}}, {center -> 6, centerIndex -> {2, 3}, 
   neighbors -> {2, 3, 1, 5, 4, 8, 9, 7}, 
   neighborsIndices -> {{1, 2}, {1, 3}, {1, 1}, {2, 2}, {2, 1}, {3, 
      2}, {3, 3}, {3, 1}}}}, {{center -> 7, centerIndex -> {3, 1}, 
   neighbors -> {6, 4, 5, 9, 8, 3, 1, 2}, 
   neighborsIndices -> {{2, 3}, {2, 1}, {2, 2}, {3, 3}, {3, 2}, {1, 
      3}, {1, 1}, {1, 2}}}, {center -> 8, centerIndex -> {3, 2}, 
   neighbors -> {4, 5, 6, 7, 9, 1, 2, 3}, 
   neighborsIndices -> {{2, 1}, {2, 2}, {2, 3}, {3, 1}, {3, 3}, {1, 
      1}, {1, 2}, {1, 3}}}, {center -> 9, centerIndex -> {3, 3}, 
   neighbors -> {5, 6, 4, 8, 7, 2, 3, 1}, 
   neighborsIndices -> {{2, 2}, {2, 3}, {2, 1}, {3, 2}, {3, 1}, {1, 
      2}, {1, 3}, {1, 1}}}}} *)

neighborIndices is a new function that takes the index of the center element and produces the indices of the neighbor elements.