% <MEDIANFILTTER> - Median filters an image in-plane
%
%   [ImOut, <SigmaOut>, <noCorrPix>, <CorrPixMask>] = 
%   medianfilter(ImIn, <kernel>, <nsigma>, <SigmaIn>, <Maskimg>, <method>)
%
%   ImOut      : Filtered matrix
%   SigmaOut   : Matrix containing the errors associated with each
%                element of ImOut
%   noCorrPix  : Gives the number of matrix elements (pixels) that 
%                have been replaced 
%   CorrPixMask: Returns a mask (either 0 or 1) containing all modified
%                matrix elements  
%   ImIn       : Matrix that has to be filtered
%   kernel     : Vector kernel = [xkernel ykernel] to select neighbours,
%                minimum: [3 3], only odd numbers allowed for xkernel and
%                ykernel
%   nsigma     : Pixels are only filtered if they are more than 
%                nsigma standard deviations away from the median/ mean 
%                within the given kernel, default nsigma = 5
%                nsigma = -1 replaces each pixel selected by Maskimg with
%                the median/ mean value of its neighbours for a given 
%                kernel
%   SigmaIn    : Matrix containing the errors associated with each
%                element of ImIn
%   Maskimg    : If Maskimg is supplied, only pixels contained in the mask
%                are filtered
%   method     : 'median' - median filter (default)
%                'mean'   - mean filter
%
%   <argument> depicts an optional argument


%==========================================================================
% 
% FUNCTION: medianfilter.m 
%           ==============
%
% $Date: 2006/12/11 12:41:23 $
% $Author: herger $
% $Revision: 1.7 $
% $Source: /import/cvs/X/PILATUS/App/lib/X_PILATUS_Matlab/medianfilter.m,v $
% $Tag: $
%
%
% <MEDIANFILTTER> - Median filters an image in-plane
% 
% Author(s):            R. Herger (RH)
% Co-author(s):         C.M. Schlepuetz (CMS)
%                       P.R. Willmott (PRW)
% Address:              Swiss Light Source (SLS)
%                       Paul Scherrer Institute
%                       CH - 5232 Villigen PSI
% Created:              2005/09/29
% 
% Change Log:
% -----------
% 
% 2005/09/29 (RH):
% - start
%
% 2005/11/09 (RH):
% - tried to implement a symmetric NaN extension -> failed (median gives
%   for every vector or matrix containing at least 1 NaN the NaN as median.
% - pixels that are more than nsigma standard deviation higher and, newly,
%   also lower the median/ mean are filtered.
%
% 2006/01/19 (RH):
% - fixed some typos.
%
% 2006/02/10 (RH & PRW): 
% - fixed median and std algorithms to make them fast 
% - (old version used even number in 3rd dimension of stack 
%   which led to very long calculating times (ca. 3 s/image)
% - New routine does this in  ca. 0.07 s! 
% 
% 2006/02/17 (PRW): 
% - corrected error in definition of CorrPixMask when nsigma set to -1.
% - Also improved the documentation. 
%
% 2006/12/11 (RH):
% - added cvs tag information for first release.

%==========================================================================
%  Main function - <medianfilter>
%                  ==============

function [ImOut, SigmaOut, noCorrPix, CorrPixMask] = ...
  medianfilter(ImIn, kernel, nsigma, SigmaIn, Maskimg, method)


%----------------------
% check input arguments

% create a vector dim = [ydim xdim] that contains the size of Im
[ydim xdim]=size(ImIn);
dim=[ydim xdim];

% are there 1 to 6 input arguments?
error(nargchk(1, 6, nargin));

% check if there are 2 ... 5 args specified
if (nargin < 6)
  method='median';                      % default method: 'median'
end;
if (nargin < 5) 
  Maskimg = ones([ydim xdim]);          % If Maskimg not specifically
end;                                    % given, fill it with ones
if (nargin < 4) 
  SigmaIn = ones([ydim xdim]);          % If Sigma matrix not specifically 
end;                                    % given, fill it with ones
if (nargin < 3) 
  nsigma = 5;                           % nsigma set to 5 as default 
end;
if (nargin < 2)
  kernel = [3 3];                       % default kernel size = [3x3]
end;

% if an empty Maskimg is supplied (i.e., []), create one filled with ones 
if (isempty(Maskimg) == 1)
  Maskimg = ones([ydim xdim]);
end;

% if an empty Sigma matrix is supplied, create one filled with ones 
if (isempty(SigmaIn) == 1)
  SigmaIn = ones([ydim xdim]);
end;

% if an empty nsigma variable is supplied, set it to 5
if (isempty(nsigma) == 1)
  nsigma = 5;
end;

% if an empty x- and y-kernel is supplied, default = 3x3
if (isempty(kernel) == 1) 
  kernel = [3 3];
end;

% check if vector kernel only has 2 elements
if (numel(kernel) ~= 2)
  error(strcat('Invalid input for ''kernel'' in function ', ...
               ' medianfilter.\nUse ''help medianfilter'' for further', ...
               ' information.\n'), ...
               '');
end;

% check if xkernel and ykernel are both odd
if (~(mod(kernel(1,1), 2) && (mod(kernel(1,2), 2))))
  error('Vector ''kernel'' has to consist of ODD numbers!');
end;

% check if xkernel and ykernel > 2
if (~((kernel(1,1) > 2) && kernel(1,2) > 2))
  error('Vector kernel has to be at least [3 3].');
end;

% select method
switch lower(method);
  case 'mean'
    meth = 1;
  case 'median'
    meth = 0;
  otherwise
    error(strcat('The argument ''method'' must either be ''median'' ',...
                 ' (default) or ''mean''!'));
end;

%----------------------
% check output argument

% prepare output matrix
ImOut = zeros([ydim xdim]);

% are 1 to 4 output arguments specified?
error(nargoutchk(1, 4, nargout));

% deal with 2 ... 4 specified output arguments
if (nargout < 4)
  CorrPixMask = zeros([ydim xdim]);
end;
if (nargout < 3)
  noCorrPix = 0;  
end;
if (nargout < 2)
  SigmaOut = zeros([ydim xdim]);
end;

%% Some debug statements
% disp(size(SigmaOut));
% disp(noCorrPix);
% disp(size(CorrPixMask));

%-------------
% medianfilter

% define some variables
n = 0;                                  % counter for stack
xkernel = kernel(1);                    % no of pix in x dim of kernel
ykernel = kernel(2);                    % no of pix in y dim of kernel
xl      = 0 - ((xkernel - 1) / 2);      % left border for circshift, etc.
xr      = 0 + ((xkernel - 1) / 2);
yu      = 0 - ((ykernel - 1) / 2);
yd      = 0 + ((ykernel - 1) / 2);

% Replicate border elements according to kernel size
ImWork = padarray(ImIn, [yd xr], 'replicate', 'both');

% create and zero some working matrices
% Stack is the stack of images used to work out the mean or median 
Stack  = zeros ([size(ImWork), (xkernel*ykernel)]); % includes central pixel 
% sdStack is the stack of images used to work out the SD. It contains only
% the shifted images and not the initial image
sdStack  = zeros ([size(ImWork), ((xkernel*ykernel)-1)]); % excludes central pixel 
MedianIm = zeros (size(ImWork));
MedianImSigma = zeros (size(ImWork));
CorrPixMask = zeros (size(ImIn));
NotCorrPixMask = zeros (size(ImIn));
CorrIm = zeros (size(ImIn));

% Fill the sdStack (3D array) with the values of the neighbouring pixels 
% along the 3rd dimension, but NOT with the real value at this position. 
% There are therefore (= xkernel*ykernel-1) images in sdStack.
for k = xl:xr
  for l = yu:yd
    if ((k==0) && (l==0))
    else
      n=n+1;
      sdStack(:,:,n) = circshift(ImWork, [l k]);
    end;
  end;
end;

Stack = sdStack;
Stack(:,:,n+1) = ImWork; % Now add the original image to sdStack to make Stack 
% Do a median/mean filter along the 3. dimension of Stack 
if (meth == 0)
  MedianIm = median(Stack, 3); % image containing all median values
else
  MedianIm = mean(Stack, 3); % image containing all mean values 
end;
% Now schnippseli back the median image to its original size 
MedianIm = MedianIm((1+yd):(size(MedianIm,1)-yd), ...
                    (1+xr):(size(MedianIm,2)-xr)); 
% Calculate the standard deviation along the 3. dim of sdStack
% This has to be done for sdStack in order to avoid producing huge 
% standard deviations if the actual pixel (the one in ImWork) is hot. 

MedianImSigma = std(sdStack, 0, 3);   % 1-D vector containing all SDs
% Now schnippseli back the std image to its original size 
MedianImSigma = MedianImSigma((1+yd):(size(MedianImSigma,1)-yd), ...
                    (1+xr):(size(MedianImSigma,2)-xr)); 

clear Stack sdStack;

% create a mask for the pixels to be corrected and count them
if (nsigma==-1) % correct whole image
    CorrPixMask = Maskimg;
else % sets those pixels which deviate by more than nsigma*SD to unity
    [CorrPixMask] = ((ImIn .* Maskimg) > ... 
      (MedianIm+nsigma*MedianImSigma)) + ...
      ((ImIn .* Maskimg) < (MedianIm-nsigma*MedianImSigma));
end;
noCorrPix = nnz (CorrPixMask);

% multiply MedianIm and errors elementwise by CorrPixMask
CorrIm = MedianIm .* CorrPixMask; % fills those pixels that deviate by 
                                  % more than nsigma*SD with the median 
                                  % values instead
CorrImSigma = MedianImSigma .* CorrPixMask;

% mask all not-corrected pixels
NotCorrPixMask = ~CorrPixMask;

% Create ImOut and SigmaOut
[ImOut] = (ImIn .* NotCorrPixMask) + CorrIm;
[SigmaOut] = (SigmaIn .* NotCorrPixMask) + CorrImSigma;

% Clear some variables
clear dim n xkernel ykernel xl xr yu yd;
clear ImWork MedianIm MedianImSigma NotCorrPixMask CorrIm CorrImSigma;

%% Some debug statements
% disp (ImIn);
% disp (ImWork);
% disp (MedianIm);
% disp (MedianImSigma);
% disp (MedianIm+nsigma*MedianImSigma);
% disp (CorrPixMask(1:100, 1:10));
% disp (noCorrPix);
% disp (NotCorrPixMask);
% disp (CorrIm);
% disp (CorrImSigma);
% disp (ImOut);
% disp (SigmaOut);

%==========================================================================
%
%---------------------------------------------------%
% emacs setup:  force text mode to get no help with %
%               indentation and force use of spaces %
%               when tabbing.                       %
% Local Variables:                                  %
% mode:text                                         %
% indent-tabs-mode:nil                              %
% End:                                              %
%---------------------------------------------------%
%
% $Log: medianfilter.m,v $
% Revision 1.7  2006/12/11 12:41:23  herger
% cvs tag information added
%
% Revision 1.6  2006/02/17 16:03:04  willmott
% changed definition of CorrPixMask for nsigma=-1 and improved documentation
%
%% Revision 1.5  2006/02/10 14:20:56  willmott
% changed median and std to improve speed
%
% Revision 1.4  2006/02/10 13:50:52  willmott
% changed median and std to improve speed
%
% Revision 1.3  2006/01/19 14:33:19  herger
% removed some typos
%
% Revision 1.2  2005/11/09 08:47:13  herger
% lower nsigma threshold implemented
%
% Revision 1.1  2005/10/24 16:34:16  herger
% matlab function: filters an image in-plane with a median or a mean filter
%
%
%================================= End of medianfilter.m ==================