function a_median = find_median_neighbours(a, n_filt)
% FIND_MEDIAN_NEIGHBOURS finds a local median for values in a matrix
%
%   A_MEDIAN = FIND_MEDIAN_NEIGHBOURS(A, N_FILT)
%
% If A is a matrix and N_FILT is a scalar then A_MEDIAN will a matrix of the
% same size as A. Each element of A_MEDIAN will be be the median of the
% square of length 2*N_FILT + 1 that is centred on the corresponding element
% in A. The exception to this is the 'outer edge' of the matrices where this
% square is not defined; in this case the elements of A and A_MEDIAN are
% identical.
%
% The algorithm is memory intensive - a 3-d array is formed by stacking
% offset versions of A on top of each other and then taking the median. To
% avoid memory usage problems for large matrices A is subdivided so that the
% calculation is done in sections.

% $Id: man133 Exp man133 $
% Copyright J. V. Mansbridge, CSIRO, Wed May  1 11:35:15 EST 2002

if (nargin ~= 2) | (nargout ~= 1)
  error('There must be 2 input arguments and 1 output argument')
end
nn = 2*n_filt + 1;
nn2 = nn*nn;
[m, n] = size(a);
if m*n < 1e6/nn2
  % disp('no split')
  m_st = 1;
  m_fi = m;
  n_st = 1;
  n_fi = n;
  if (nn >= m) | (nn >= n)
    error('n_filt is too big for this algorithm and the given matrix')
  end
else
  % disp('split')
  m_sub = ceil(m/nn);
  m_st = [1 (m_sub*(1:(nn-1)) - n_filt)];
  m_fi = [(m_sub*(1:(nn-1)) + n_filt) m];
  n_sub = ceil(n/nn);
  n_st = [1 (n_sub*(1:(nn-1)) - n_filt)];
  n_fi = [(n_sub*(1:(nn-1)) + n_filt) n];
  if any(m_st < 1) | any(n_st < 1) | any(m_fi > m) | any(n_fi > n)
    error('n_filt is too big for this algorithm and the given matrix')
  end
end
a_median = a;
for ii = 1:length(m_st)
  for jj = 1:length(n_st)
    m_big = m_fi(ii) - m_st(ii) + 1 - 2*n_filt;
    n_big = n_fi(jj) - n_st(jj) + 1 - 2*n_filt;
    if (m_big < 1) | (n_big < 1)
      error('n_filt is too big')
    end
    a_sub = a(m_st(ii):m_fi(ii), n_st(jj):n_fi(jj));
    a_big = zeros(m_big, n_big, nn2);
    kk = 0;
    for ik = 0:2*n_filt
      for jk = 0:2*n_filt
	kk = kk + 1;
	a_big(:, :, kk) = a_sub((1+ik):(m_big+ik), (1+jk):(n_big+jk));
      end
    end
    a_median((m_st(ii) + n_filt):(m_fi(ii) - n_filt), ...
	     (n_st(jj) + n_filt):(n_fi(jj) - n_filt)) = median(a_big, 3);
    
  end
end
