source: prextra/neighbh.m

function J = neighbh(a,nbtype)
% J = neighbh(sz,nbtype)
% J = neighbh(a,nbtype)
% 
% Compute the indices of the pixels that are in the neighborhood of a
% pixel. This is done for all pixels in a whole image and returned in
% the index array J. The pixels are enumerated in this standard order
% (for an 3x5 image):
%     1   4   7  10  13
%     2   5   8  11  14
%     3   6   9  12  15
% When there is no neighbor (i.e. the pixel is on the boundary), index 0
% is returned.
% Then J becomes for a '4NN' neighborhood:
%     0  0  2  4      % pixel 1 have neighbors 2 and 4...
%     0  1  3  5      % pixel 3 has 3 neighbors..
%     0  2  0  6
%     1  0  5  7
%     2  4  6  8
%     and so on ...
%
% We can define neighborhood types nbtype:
%   '4NN'   4-nearest neighborhood
%   '8NN'   8-nearest neighborhood
%   '9NN'   9-nearest neighborhood
%  '25NN'  25-nearest neighborhood
% It is also possible to give your own definition of the neighborhood.
% You should give the indices w.r.t. the origin, for example:
%   nbtype = [-3 0;
%             +3 0;
%              0 -1];
%

if isdataset(a)
        if ~chckcoord(a)
                error('Does not look like an image dataset');
        end
        sz = getobjsize(a);
else
        if ~(size(a)==[1 2])
                error('I want to have the size of the image');
        end
        sz = a;
end
% The total number of pixels in the image:
N = prod(sz);
% The definitions of the predefined neighborhoods:
if ischar(nbtype)
        switch nbtype
        case '4NN'
                % left    up     down    right
                pos = [-1 0; 0 -1; 0  1; 1 0];
        case '8NN'
                pos = [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1];
        case '9NN'
                pos = [-1 -1; -1 0; -1 1; 0 -1; 0 0; 0 1; 1 -1; 1 0; 1 1];
        case '25NN'
                pos = [-2 -2; -2 -1; -2  0; -2  1; -2  2; -1 -2; -1 -1;...
                         -1  0; -1  1; -1  2; 0 -2; 0 -1; 0  0; 0  1; 0  2; ...
                          1 -2; 1 -1; 1  0; 1  1; 1  2; 2 -2; 2 -1; 2  0;...
                          2  1; 2  2];
        case '4star'
                %      left   up   down right
                pos = [-2 0; 0 -2; 0  2; 2 0];
        case '4cross'
                %      upleft  upright downleft downright
                pos = [-2 -2;   2 -2;   -2  2;   2 2];
        otherwise
                error('This neighborhood is not defined');
        end
else
        if size(nbtype,2)~=2
                error('The neighborhood definition should be a nx2 matrix.');
        end
        pos = nbtype;
end
% convert these neighborhood positions into changes in indices:
dx = sz(1)*pos(:,1)' + pos(:,2)';
% in horizontal limits (lower and upper):
ckhorz = -pos(:,1)';
I = find(ckhorz<0);
ckhorz(I)=ckhorz(I)+sz(2)+1;
% and in vertical limits (lower and upper):
ckvert = -pos(:,2)';
I = find(ckvert<0);
ckvert(I)=ckvert(I)+sz(1)+1;

% Setup the general parameters
D = length(dx);

% Define the pixels, and compute their coordinates:
I = (1:N)';
hcoord = floor((I-1)/sz(1))+1;
vcoord = mod(I-1,sz(1))+1;

% here the results will be stored:
J = zeros(N,D);

% and here we go:
for i=1:D
        % Here the main work of constructing the J matrix is performed,
        % it is actually just these two lines. The trick is to shift the I
        % vector up or down and store it alongside the already created matrix
        % J.

        if dx(i)>0   % look at pixels in the + direction

                J(:,i) = [I((1+dx(i)):end); zeros(dx(i),1)];

        else         % look at pixels in the - direction

                J(:,i) = [zeros(-dx(i),1); I(1:end+dx(i))];

        end
        % Unfortunately, most of the work is actually in the checking of the
        % bounds:
        % Check for the horizontal bounds:
        % look at the pixels left:
        if pos(i,1)<0
                if ckhorz(i)
                        excp = find(hcoord<=ckhorz(i));
                        J(excp,i) = 0;
                end
        else  % and pixels right
                if ckhorz(i)
                        excp = find(hcoord>=ckhorz(i));
                        J(excp,i) = 0;
                end
        end
        % Check for the vertical bounds:
        % look at the pixels *above*:
        if pos(i,2)<0
                if ckvert(i) 
                        excp = find(vcoord<=ckvert(i));
                        J(excp,i) = 0;
                end
        else   % otherwise the pixels below:
                if ckvert(i)
                        excp = find(vcoord>=ckvert(i));
                        J(excp,i) = 0;
                end
        end
end

return


function out = chckcoord(a)
% out = chckcoord(a)
%
% check if dataset a contains pixels regularly sampled from an image
% It looks a the objsize and maybe the X- and Y-coordinates.

out = 1;

sz = getobjsize(a);

% first check if we even have some 2D objsize
if length(sz)==1
        out = 0;
        return
end

featlab = getfeatlab(a);
% now check the possible X coordinate:
xcoord = strmatch('X',featlab);
if ~isempty(xcoord)
        cx = +a(:,xcoord);
        dx = diff(cx(1:sz(1)));
        if std(dx)~=0
                out = 0;
                return
        end
end
% now check the possible Y coordinate:
ycoord = strmatch('Y',featlab);
if ~isempty(ycoord)
        cy = +a(:,ycoord);
        dy = diff(cy(1:sz(1):end));
        if std(dy)~=0
                out = 0;
                return
        end
end

return