% <HISTMASK> - Produces a mask based on two threshold values
%
%
%   [MaskImg, <noPix>, <noLowHigh>] = histmask(Im, thresh, <MaskimgIn>)
%
%   MaskImg  : MaskImg based on Im. MaskImg is
%              1 for all pixels between or equal low and high,
%              0 otherwise
%   noPix    : Returns the number of pixels within the boundaries
%   noLowHigh: Vector containing no of pixels lower and higher
%              than thresh. noLowHigh = [Low High]
%
%   Im       : Matrix or stack of matrices to mask
%   thresh   : Vector containing both the lower (low) and higher (high)
%              threshold, thresh = [low high]. Note: apixel = low or 
%              apixel = high will be masked as 1.  
%   MaskImgIn: If MaskimgIn is supplied, only pixels contained in the mask
%              are processed, e.g. if you want to filter a region 
%              of interest.
%              Note that Im and MaskimgIn have to be of the same size in
%              xdim and ydim. zdim must either be the same or 1
%              to avoid confusion.
%
%   <argument> depicts an optional argument

%==========================================================================
% 
% FUNCTION: histmask.m 
%           ==========
%
% $Date: 2006/12/11 12:39:13 $
% $Author: pauli_s $
% $Revision: 1.5 $
% $Source: /import/cvs/X/PILATUS/App/lib/X_PILATUS_Matlab/histmask.m,v $
% $Tag: $
%
%
% <HISTMASK> - Produces a mask based on two threshold values
% 
% Author(s):            R. Herger (RH)
% Co-author(s):         C.M. Schlepuetz (CS), D. Martoccia (DM)
% Address:              Swiss Light Source (SLS)
%                       Paul Scherrer Institute
%                       CH - 5232 Villigen PSI
% Created:              2004/09/16
% 
% Change Log:
% -----------
% 
% 2005/11/02 (RH)
% - created
%
% 2006/10/20 (SAP)
% - corrected orange errors
% - the function works also with a stack of images


%==========================================================================
%  Main function - histmask
%                  ========

function [MaskImg, noPix, noLowHigh] = histmask(Im, thresh, MaskImgIn)

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

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

% are there 2 to 3 input arguments?
error(nargchk(2, 3, nargin));

if (nargin < 3) 
  MaskImgIn = ones([ydim xdim zdim]);    % Maskimg matrix filled with ones
end;

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

% if empty MaskImgIn is supplied, create a MaskimgIn with 1
if (isempty(MaskImgIn) == 1)
  MaskImgIn = ones([ydim xdim zdim]);
end;

% If the dimenstion of MaskImgIn in z-direction is just 1 and the dimension
% of Im in z-direction is not, copy the entries to make the arrays
% equalsized.
if (size(MaskImgIn,3) ~= zdim) && (size(MaskImgIn,3) == 1)
    MaskImgIn = padarray(MaskImgIn,[0 0 zdim-1],'replicate','post');
end;

% check if dimensions of Im and of MaskImgIn are the same
if (~isequal(size (MaskImgIn), [ydim xdim zdim]) && ...
        ~isequal(size (MaskImgIn), [ydim xdim]))
  error('xdim, ydim and zdim of Im and MaskImgIn must agree!');
end;


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

% prepare output matrix default to unity
MaskImg = ones([ydim xdim zdim]);

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

%---------
% histmask

% apply MaskimgIn to Im
ImWork = Im .* MaskImgIn;


% create indices for pixels below and above the thresholds
indicesLo = (ImWork < thresh(1)); % vector list of low bad pix posns
indicesHi = (ImWork > thresh(2)); % vector list of hi bad pix posns

% count pixels below and above thresh and create vector noLowHigh
[noLowHigh(:,1)] = arrayfun(@(x)nnz(indicesLo(:,:,x)),1:zdim);
[noLowHigh(:,2)] = arrayfun(@(x)nnz(indicesHi(:,:,x)),1:zdim);

% set Maskimg below and above thresh to zero
MaskImg(indicesLo) = 0;
MaskImg(indicesHi) = 0;

% find nonzero elements
% indicesMask = find (ImWork ~= 0);

%  count nonzero elements
[noPix] = arrayfun(@(x)nnz(MaskImg(:,:,x)),1:zdim);

% Clear some variables
clear dim indicesLo indicesHi ImWork


%==========================================================================
%
%---------------------------------------------------%
% 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: histmask.m,v $
% Revision 1.5  2006/12/11 12:39:13  pauli_s
% added tag
%
% Revision 1.4  2006/11/30 13:57:56  pauli_s
% Remove unneeded files
%
% Revision 1.3  2006/02/15 16:17:23  willmott
% removed bug
%
% Revision 1.2  2005/11/04 09:13:21  herger
% removed typo
%
% Revision 1.1  2005/11/02 16:07:22  herger
% matlab function: produces a mask based on two threshold values
%
%
%================================= End of histmask.m ======================