function ngb = get4ngb(rows,cols,x,y)
  
  % function ngb = get4ngb(rows,cols,x,y)
  % This function returns the cooridinates of the 8-neighbours
  % of an element in a matrix.  The values are returned in the
  % list ngb.  If the neighbour does not exist,- that is (x,y)
  % corredsponds to an edge or corner of the array,
  % the list of neighbours contains only those coordinates corresponding
  % to real neighbours.
  % Copyright Lars Aurdal/FFIE.
  
  % Handle left edge.
  
  % Handle upper left corner.
  
  if ((x ==  1) & (y == 1))
    ngb(1:2,1:1) = [2 1]';
    ngb(1:2,2:2) = [1 2]';
    
    % Handle lower left corner.
    
  elseif ((x == rows) & (y == 1))
    ngb(1:2,1:1) = [rows 2]';
    ngb(1:2,2:2) = [(rows - 1) 1]';
    
    % Handle left edge in general.
    
  elseif (y == 1)
    ngb(1:2,1:1) = [(x+1) 1]';
    ngb(1:2,2:2) = [x 2]';
    ngb(1:2,3:3) = [(x-1) 1]';
    
    % Handle right edge.
    
    % Handle upper right corner.
    
  elseif ((x == 1) & (y == cols))
    ngb(1:2,1:1) = [1 (cols-1)]';
    ngb(1:2,2:2) = [2 cols]';
    
    % Handle lower right corner.
    
  elseif ((x == rows) & (y == cols))
    ngb(1:2,1:1) = [rows (cols-1)]';
    ngb(1:2,2:2) = [(rows-1) cols]';
    
    % Handle right edge in general.
    
  elseif (y == cols)
    ngb(1:2,1:1) = [(x+1) cols]';
    ngb(1:2,2:2) = [x (cols-1)]';
    ngb(1:2,3:3) = [(x-1) cols]';
    
    % Handle top line.
    
  elseif (x == 1)
    ngb(1:2,1:1) = [1 (y-1)]';
    ngb(1:2,2:2) = [2 y]';
    ngb(1:2,3:3) = [1 (y+1)]';
    
    % Handle bottom line. 
    
  elseif (x == rows)
    ngb(1:2,1:1) = [rows (y-1)]';
    ngb(1:2,2:2) = [(rows-1) y]';
    ngb(1:2,3:3) = [rows (y+1)]';
    
    % Handle general case
    
  else
    ngb(1:2,1:1) = [(x-1) y]';
    ngb(1:2,2:2) = [x (y-1)]';
    ngb(1:2,3:3) = [(x+1) y]';
    ngb(1:2,4:4) = [x (y+1)]';
    
  end