How to find Moore neighbourhood in a 2D cellular automata?

I am trying to find the list of neighbours inside a Moore neighbourhood of variable length in a 2D cellular automata (will run different conditions on them later on). I have attached a link of an image for a Moore neighbourhood with variable r from Wolfram MathWorld, in case anyone is interested. In an NXN grid, my code is as below. This is ok for r = 1; but for r>1, it keeps missing some values. For example, for r = 2, it misses (i-1,j-2);(i+1,j-2);(i+2,j-1);(i+2,j+1);(i+1,j+2);(i-1,j+2);(i-2,j+1);(i-2,j-1).

for i = 1:N
  for j = 1:N
     for r = 1:R
       neighbour_ij = [(i-r,j-r);(i,j-r);(i+r,j-r); (i+r,j);(i+r,j+r);(i,j+r);(i-r,j+r);(i-r,j)];
  end
end

Answer 1

[Bee]

In case someone is looking for this, I am sharing the way I solved it - there are more efficient ways to do this, which I don't know, but mine at least works :)

% Take a grid, G = NXN
% define a neighbourhood matrix
moore = cell(size(Z));
% for example, r = 1;  the neighbourhood size  = (2*r+1)^2

for i = 1:N
   for j = 1:N
   x = [i-r:i+r]; %row subscript 
   x(x<=0)=[]; x(x>N)=[]; %discarding values that are out of range, for absorbing boundary
      %for toroidal neighbourhood, discard the previous line which is for absorbing boundary
      %x(x<=0) = N+ x(x<=0);
      %x(x>N) = -N+ x(x>N);
      %end of toroidal neighbourhood
   y = [j-r:j+r]; %column subscript
   y(y<=0) = []; y(y>N) = []; % absorbing boundary, i.e. edge nodes have fewer neighbours, corner 
                              % node has only 3 neighbours etc.
      %for toroidal neighbourhood, discard the previous line which is for absorbing boundary
      %x(y<=0) = N+ y(y<=0);
      %x(y>N) = -N+ y(y>N);
      %end of toroidal neighbourhood

     for k = 1:length(x)
        for l = 1:length(y)
            moore{i,j} = [moore{i,j};(Z(x(k),y(l)))]; 
        end
     end
     moore{i,j} = unique(moore{i,j});
     end
  end

Answer 2

[KSSV]

Read about knnsearch

[Bee]

I couldn't relate them, knnsearch finds nearest neighbours in a matrix, X comparing it with another matrix. Y - that's what I understood from reading the doc, whereas I am looking for nearest neighbour of each element in vector X.