function [Mout, badi] = medfilt2d(M, vert_filtlength, horiz_filtlength, tol);
% [Mf, badi] = medfilt2d(M, vert_filtlength, horiz_filtlength, tol);
%
% M is a matrix (eg.  50xN timeseries of raw adcp data)
% vert_filtlength and horiz_filtlength must be odd
% values which are tol away from the median are replaced by the median
%
% NOTE: horizontal loop is executed first
%
% given a big matrix, replace the middle value (middle depth and middle 
%   profile) with the median from that comatrix
% eg: 
% [Mf, badi] = medfilt2d(amp, 11, 5, 10);


%looking at it as if it were (50xN0 timeseries of raw adcp data
%choose how many profiles wide and how many bins deep to make the window

%this is just a hack of medfilt

if (nargin ~= 4)
  help medfilt2d
  error('requires 4 arguments')
end

%require vert_filtlength and horiz_filtlength to be odd
if rem(vert_filtlength,2) ~= 1
  error('vertical filter length must be odd-valued')
end
if rem(horiz_filtlength,2) ~= 1
  error('horizontal filter length must be odd-valued')
end

if (  (vert_filtlength == 1) & ( horiz_filtlength == 1) )
  error('why bother!! ')
end


% eg (8x13)
[nbins, nprofs] = size(M);

%pad the left and right with the correct number of profiles
% eg: vert_filenum = 5; horiz_filtlength = 3;
% this means  half_vfiltlength = 2; half_hfiltlength = 1;

half_vfiltlength = fix(vert_filtlength/2);%pad top and bottom by this many bins
half_hfiltlength = fix(horiz_filtlength/2);%pad sides by this many profiles


% padlr would be [M(:,1) M(:,1)]  in this case (i.e. two profiles = padding)
padl = (ones(half_hfiltlength, 1) * M(:,1)')';
padr = (ones(half_hfiltlength, 1) * M(:,end)')';
paddedM = [padl M padr];


%used to be for a row vector (i.e. column')  
% so this would be one vertical column (transposed) repeated horiz_filt times
%w = NaN * ones(fl,ncol+fl-1);
%      ones(vert_filtlength, nbins + vert_filtlength - 1)
%      all stacked up as rows now, with profiles going into filtered
%      group stacked up in the vertical

%ick.  loop through the profiles 
Mmed = NaN *ones(size(M));
for profi = 1:(nprofs)  %loop through each profile
  grab_profi = profi + [1:horiz_filtlength] - 1; % (1,2,3), (2,3,4), ...
  
  w = NaN * ones(vert_filtlength * horiz_filtlength,...
	  nbins + vert_filtlength -1);
  % for ii = 1:fl ;   w(ii,ii:ii+ncol-1) = f; end
  % for ii = 1:vert_filtlength  ;   %  w(ii, ii+nbins-1) = profile; end
   
  for vcounti = 1:vert_filtlength
	%indices into w:
	wii = (vcounti - 1)*horiz_filtlength + [1:horiz_filtlength];
	w(wii,vcounti + [1:nbins] - 1) = paddedM(:,grab_profi)';
  end
  %from medfilt:
  %m = median_nan(w(:,half_fl+1:half_fl+ncol));
  %here this goes in the next column of Mout
  [j_mean, j_std, j_n, j_med]  = ...
	  mstdgap(w(:,half_vfiltlength+1:half_vfiltlength+nbins));
  Mmed(:,profi) = j_med';
  
end

Mout = M;

badi = find(abs(Mout-Mmed)>tol);
if (~isempty(badi))
  Mout(badi)=Mmed(badi);
end
