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

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).

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

partitionMap[f_, m_, r_] := DeveloperPartitionMap[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}}}} *)


DeveloperPartitionMap 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 DeveloperPartitionMap 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}

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.

• Very nice! A couple of questions: what part of this solution takes care of the periodic boundaries, and what does DeveloperPartitionMap does? (and, it seems that it has been "superseded" by BlockMap', are they equivalent?) Thanks! Commented Jul 3, 2020 at 23:54
• It seems that the expected number of neighbors in a Moore/von Neumann radius of 2 (or more) does not match the solution here. Maybe they are also swapped (i.e. for a Moore neighborhood of radius 2 one expects 24 neighbors, but we get only 8). The "Moore" option seems to give the Neumann neighborhoods, but undercounted (i.e. a von Neumann of radius 2 should give 12 neighbors, but it only gives 10 -when using "Moore", as I mention-). Commented Jul 4, 2020 at 0:57
• @TumbiSapichu I've corrected the problems now, please have another look. I also made a comment about PartitionMap, that it takes care of the periodic boundaries using arguments documented under the Partition function. I'm not sure how the indexing is supposed to be done, do you mean like the index that that particular neighbor has in the full matrix? Or its index in the neighborhood? Commented Jul 4, 2020 at 7:28
• @TumbiSapichu Have a look at my new update. Commented Jul 4, 2020 at 20:34
• @TumbiSapichu They are not guaranteed to continue to be supported, although I would say they are likely to be. You are right that partitionMapIndexed does not rely on that function, and you may also note that if you replace MapIndexed in that function by Map then you get the exact equivalent of DeveloperPartitionMap for matrices. Commented Jul 4, 2020 at 20:47