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
