% FFTEST - Tests the Flatfield Data and generates the Flatfield 
% correction image with the standard deviation of each pixel
%
%   [Im , <sigma>, <info>] = fftest(ImStack, <fileOut>, <printer>)
%
%   Im     : Flatfield correction image, normalized around 1 to correct
%            a pixel's' different quantum efficiency for incoming photons. 
%   sigma  : Standard deviation of each pixel of Im (ideally, this reflects
%            the counting statistics of each pixel, normalized to the mean
%            of Im).
%   info   : struct-variable with additionanal information of the analysis
%            of the flatfield: 
%               info.mean : Mean value of all flatfields
%               info.low  : lower count-acceptance tolerance
%               info.up   : upper count-acceptance tolerance
%               info.tol  : tolerance percentage of outliers
%               info.bad  : number of bad pixels
% 
%   ImStack : Flatfield data, generated by fillstack
%   fileOut : optional output-"mat"-file, with the fields
%                 ff_dat.ff_mask = Im;
%                 ff_dat.ff_err  = sigma;
%                 ff_dat.ff_info = info;
%             default is ff_corr.mat
%   printer : optional argument which specifies the printer. default is
%             WBBA_005_1
% 
%    <argument> depicts an optional argument

%==========================================================================
% Further info for function FFTEST
% --------------------------------
% Function subroutine:
% -------------------
%
% HISTSUMMASK - gets the histmask from the summed histmasks of each and
% every ff image taken. 
%
% POISSON - calculates a theoretical poisson distribution and delivers the 
%           value at position x
%

%==========================================================================
% 
% FUNCTION: fftest.m 
%           ==========
%
% $Date: 2006/12/11 12:39:40 $
% $Author: pauli_s $
% $Revision: 1.2 $
% $Source: /import/cvs/X/PILATUS/X04SA/App/lib/X_PILATUS_X04SA_Matlab/fftest.m,v $
% $Tag: $
%
% 
% Author(s):            S. Pauli (SP)
% Co - Authors          Roger Herger (RH), Christian Schlepuetz (CS)
% Address:              Swiss Light Source (SLS)
%                       Paul Scherrer Institute
%                       CH - 5232 Villigen PSI
% Created:              2006/02/20
% 
% Change Log:
% -----------
% 
% 2006/02/20 (SP)
% - created
%
% 2006/02/22 (SP)
% - changed Input to just ImStack, instead of the same variables, as
%   fillstack has
%
% 2006/03/09 (SP)
% - added printing and saving dialog and set the plots on 2 figures
%  
% 2006/04/05 (SP)
% - optimized input and presentation, axis of first histogram are now 
%   changeable, output will be saved in a file
%
% 2006/05/06 (SP)
% - Saved file is a mat-file, mask and maskerror as struct, added printer
% input in input arguments
%
% 2006/12/07 (SP)
% - Made the graphical output nicer
% - Added a theoretical poissonian into the histogram plot
% - added output-variable info

%==========================================================================
%  
%% Main function - fftest
%                  ======
function [Im sigma info]=fftest(ImStack,fileOut,printer)

%% check input arguments

% is there 1 input arguments?
error(nargchk(1, 3, nargin));

if (nargin < 3)
    printer = 'WBBA_005_1' ; 
   %printer = 'WSLA_X04_2' ;
   %printer = 'WLGA_227_1' ;
   if (nargin < 2)
       fileOut='ff_corr.mat';
   end;
end;

% check if ImStack has 3 dimensions
if (ndims(ImStack) ~= 3)
  error(strcat('Invalid input for ''Stack'' in function ', ...
               ' immerge.\nUse ''help immerge'' for further', ...
               ' information.\n'), ...
               '')
end;

% Some variables read out of the input variables
[ydim xdim zdim] = size(ImStack);

%% check output argument 

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

%% Calculating and plotting of a histogram 
%    with counting statistic of each pixel of the stack

fprintf(1,'Preparing histogram...');

set(0,'DefaultFigureWindowStyle','docked');
fig1 = figure;

% ---
% preparing of the stack and making a histogram
shstack = reshape(ImStack,1,[])'; % vector with cts per pixel
peak = mean(shstack); % mean value of shstack
xmax = 3*round(peak); % the maximal value of histogram bins
bin = 0:1:xmax;
[hdat,bdat]=hist(shstack,bin);

%----
% find the position of the histogram's maximal value
m=max(hdat(2:size(hdat,2)-1)); 
peak = find(hdat==m);
sigp = sqrt(peak);

axes('units','normalized','position',[0.06 0.72 0.38 0.22]);
    bar(bdat,hdat);
    grid on;    hold on;
    xlim([-2 xmax]);
    xlabel('#counts');
    ylabel('#pixels');
    plot(1:1:xmax , ...
        length(shstack)*poisson(1:1:xmax,find(hdat==max(hdat))),'r');
    title({'Average counting-rate of each pixel';...
        [' over ',num2str(zdim),' images']});
    ylim([0 1.1*max( ...
        length(shstack)*poisson(1:1:xmax,find(hdat==max(hdat))) )]);

clear shstack bdat hdat;
%% Input dialog for changing the axis

fprintf(1,'done\n');
axReply = input(['Do you want to change the ',...
    'limits of the x- and y-axes? y/n [n]: '],'s');
if isempty(axReply)
   axReply = 'n';
end;

ylimits = get(gca,'ylim');
  
if (axReply=='y')

    reply = 'n';

    axmin = -2;
    axmax = xmax;
    aymax = ylimits(2);

    while (reply ~='y')

        rpAxmin = ...
            input(['Set lower x - limit  [',int2str(axmin),']: ']);
        if ~isempty(rpAxmin)
            axmin = rpAxmin;
        end;

        rpAxmax = ...
            input(['Set upper x - limit  [',int2str(axmax),']: ']);
        if ~isempty(rpAxmax)
            axmax = rpAxmax;
        end;
        while (axmax < axmin)
            rpAxmax = ...
                input(['Set UPPER!!! x - limit  [',int2str(axmax),']: ']);
            if ~isempty(rpAxmax)
                axmax = rpAxmax;
            end;
        end;

        rpAymax = ...
            input(['Set upper y - limit  [',int2str(aymax),']: ']);
        if ~isempty(rpAymax)
            aymax = rpAymax;
        end;

        axis([axmin axmax 0 aymax]);

        reply = input('Are the axis OK? y/n [y]: ', 's');
        if isempty(reply)
            reply = 'y';
        end;

    end;
end;

info.mean = peak;
%% Input dialog for the limits of the histogram 

reply = 'n';

loLim = round((peak-6*sigp)/10)*10;
hiLim = round((peak+6*sigp)/10)*10;


while (reply ~='y')

    rpLoLim = input(...
        ['Set lowest acceptable counting level [',int2str(loLim),']: ']);
    if ~isempty(rpLoLim)
        loLim = rpLoLim;
    end;

    line([loLim,loLim],[0,ylimits(2)],'Color','b');

    rpHiLim = input(...
        ['Set highest acceptable counting level [',int2str(hiLim),']: ']);
    if ~isempty(rpHiLim)
        hiLim = rpHiLim;
    end;

    line([hiLim,hiLim],[0,ylimits(2)],'Color','b');

    reply = input('Are these values OK? y/n [y]: ', 's');
    if isempty(reply)
        reply = 'y';
    elseif (reply ~='y')
        line([loLim,loLim],[0,ylimits(2)],'Color',[0.5 0.5 0.5]);
        line([hiLim,hiLim],[0,ylimits(2)],'Color',[0.5 0.5 0.5]);
    end;


end;

info.low = loLim;
info.up  = hiLim;

clear peak;
%% Calculation of the Tolerance and plot the # of bad pix

fprintf(1,'Calculating Tolerance Graph...\n');

% Calculation of the Tail

[maskimgout, summedmask, nbp] = ...
    histsummask(ImStack, [loLim hiLim], 1);

tol=zdim;
j=1;

pixct = 0;
n     = 0;

while (tol>0)

    [maskimgout, summedmask, nobadpix] = ...
        histsummask(ImStack, [loLim hiLim], tol);
    pixct(j) = nobadpix ;
    n(j) = zdim-tol ;
    fprintf(1,strcat(int2str(zdim-tol),', '));
    tol=tol-1;
    if (nobadpix ~=0)
        if (tol < zdim-5) && (nobadpix < 2*nbp) && ...
                ( ( (oldnobadpix-nobadpix) /nobadpix) <0.25) && (tol~=0)
            if  (tol>20)
                if (tol<zdim/2)
                    tol=1;
                else
                    tol=tol-9;
                end;
            else
                tol=1;
            end;
        end;
    end;

    j=j+1;
    oldnobadpix=nobadpix;

    if (tol==0)
        [maskimgout, summedmask, nobadpix] = ...
            histsummask(ImStack, [loLim hiLim], 0);
        pixct(j) = nobadpix ;
        n(j) = zdim-tol ;
        fprintf(1,strcat(int2str(zdim-tol),'\n'));

    end;

end;

axes('units','normalized','position',[0.56 0.72 0.38 0.22]);
plot((zdim-n)/zdim*100,pixct);
      view(180,-90);
      grid on;
      title('Incidence of bad pixels vs. tolerance');
      ylimits = get(gca,'ylim');
      xlimits = get(gca,'xlim');
      yl = (ylimits(2)-ylimits(1))/20;
      text(xlimits(1),ylimits(2)-yl,{['total no. of pixels: ', ...
          int2str(xdim*ydim)];'Definition of tolerance:'; ...
          'e.g. 80% means 4 out of 5 times';...
          'the pixel lies within the ';...
          'acceptable counting rates'},...
          'HorizontalAlignment','right', ...
          'VerticalAlignment','top');       
      xlabel('Tolerance');
      ylabel('Bad Pixels');
%% Input dialog for the limit of tolerance

reply = 'n';
tolLim = 98;

while (reply ~='y')
  
   rpTolLim = ...
       input(['Set Tolerance limit in % [',int2str(tolLim),'%]: ']);
   if ~isempty(rpTolLim)
      tolLim = rpTolLim;
   end;
   
   line([tolLim,tolLim],ylimits,'Color','r');
    
   reply = input('Is that value OK? y/n [y]: ', 's');
   if isempty(reply)
       reply = 'y';
   else   
      line([tolLim,tolLim],ylimits,'Color',[0.5 0.5 0.5]);
   end;

end;

clear pixct;

info.tol = tolLim;
%% Calculating the mask with the user-given tolerance

fprintf(1,'Calculating the mask...');

[MaskImgOut, summedmask, nobadpix] = ...
    histsummask(ImStack, [loLim hiLim], round(tolLim/100*zdim));
info.bad = nobadpix;

fprintf(1,'done\nCalculating the corrected flatfield and error...');
[Im,sigma] = makeffcorr(ImStack,MaskImgOut);

clear summedmask;

ff_dat.ff_mask = Im;
ff_dat.ff_err  = sigma;
ff_dat.ff_info = info;
save(fileOut,'-struct','ff_dat');

%% Calculating the std.deviation of one pixel and plot the deviationfield

fprintf(1,'done\nCalculating z-direction standard-deviation...');

timeErr=std(ImStack,0,3);
sCop=mean2(timeErr);

ax1 = axes('units','normalized','position',[0.06 .39 0.6 .22]);
cTol=round(sqrt(sCop));
image(timeErr,'CDataMapping','scaled');
   set(gca,'clim',[sCop-3*cTol sCop+3*cTol]);
   axis([0 xdim 0 ydim]);
   colormap(pilatus)
   axis off;
   title('Standard Deviation (SD) of each Pixel');
   axes('units','normalized','position',[0.68 .39 0.04 .22]);
   handle = gca;
   colorbar(handle,'peer',ax1);

%% Find the worst pixel and a mean pixel and plot all their values

maxErr=max(max(timeErr.*MaskImgOut));
[indMaxErrX,indMaxErrY] = find(timeErr == maxErr);
[indMeanErrX,indMeanErrY] = ...
       find((timeErr > sCop-0.0001) & (timeErr < sCop+0.0001) );

j=-3;
while isempty(indMeanErrX) 
    [indMeanErrX,indMeanErrY] = ...
       find((timeErr.*MaskImgOut > sCop-10^(j)) & ...
       (timeErr.*MaskImgOut < sCop+10^(j)) );
    j=j+1;
end;

axes('units','normalized','position',[0.06 0.06 0.38 0.22]);
   plot(squeeze(ImStack(indMaxErrX,indMaxErrY,:)),'Color','r');
   xlim([1 zdim]);
   ax1 = gca;
   set(ax1,'XColor','r','YColor','r')

   ax2 = axes('Position',get(ax1,'Position'),'XAxisLocation','top',...
       'YAxisLocation','right','Color','none', ...
       'XColor','b','YColor','b');
   grid on;
   hold on;

   plot(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:)),...
       'Color','b','Parent',ax2);
   xlim([1 zdim]);
   title('The worst pixel and a typical pixel');
%% Plotting a histogram over the std.deviation of every pixel

edges=0:1:3*sCop+1;
[hisTimeErr binTimeErr]=hist(reshape(timeErr,[],1),edges);

axes('units','normalized','position',[0.78 0.39 0.16 0.22]);
   bar(binTimeErr,hisTimeErr);
   xlim([-1 3*sCop+1]);
   title({'Histogram of the';'SD error image (left)'});
   set(gca,'YAxisLocation','right');

clear timeErr;
%% Plots a histogram of the worst and a typical pixel

loMeanAx=fix(min(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:))));
hiMeanAx=round(max(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:))));
binMeanAx=round((hiMeanAx-loMeanAx)/50);
hisMeanTick=loMeanAx:binMeanAx:hiMeanAx;
[hMeanDat nMeanBin]=  ...
    hist(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:)),hisMeanTick);

axes('units','normalized','position',[0.56 0.06 0.16 0.22]);
   bar(nMeanBin,hMeanDat);
   axis([loMeanAx-1 hiMeanAx+1 0 max(hMeanDat)+1]);
   title('Typical Pixel');
   xlimits = get(gca,'xlim');
   set(gca,'XTick',...
       [xlimits(1) round((xlimits(1)+xlimits(2))/2) xlimits(2)]);

loMaxAx=fix(min(squeeze(ImStack(indMaxErrX(1),indMaxErrY(1),:))));
hiMaxAx=round(max(squeeze(ImStack(indMaxErrX(1),indMaxErrY(1),:))));
binMaxAx=round((hiMaxAx-loMaxAx)/50);
hisMaxTick=loMaxAx:binMaxAx:hiMaxAx;
[hMaxDat nMaxBin]=...
    hist(squeeze(ImStack(indMaxErrX(1),indMaxErrY(1),:)),hisMaxTick);

axes('units','normalized','position',[0.78 0.06 0.16 0.22]);
   bar(nMaxBin,hMaxDat);
   axis([loMaxAx-1 hiMaxAx+1 0 max(hMaxDat)+1]);
   title('Worst Pixel');
   xlimits = get(gca,'xlim');
   set(gca,'XTick',...
       [xlimits(1) round((xlimits(2)+xlimits(1))/2) xlimits(2)]);

%% Printing and saving dialog

reply = 'n';
ansPrint='';
printerNam = printer;
PSnam      = 'stats.eps';

while (reply ~='y')

    ansPrint = input('\nDo you want to print the Figure? [y] :','s');
    if isempty(ansPrint)
        ansPrint = 'y';
    end;

    if (ansPrint=='y')
        namPrint = input( ...
            ['What is the name of the printer? [',printerNam,'] : '],'s');
        if isempty(namPrint)
            namPrint=printerNam;
        end;
        namPrint2=strcat('-P',namPrint);
    end;

    ansPS = input( ...
        'Do you want to save the figure as an eps-file? [y] :','s');
    if isempty(ansPS)
        ansPS = 'y';
    end;

    if (ansPS=='y')
        namPS = input( ...
            ['How should the eps-file be called? [',PSnam,'] : '],'s');
        if isempty(namPS)
            namPS=PSnam;
        end;
    end;

    reply = input('Are all inputs OK? y/n [y]: ', 's');
    if isempty(reply)
        reply = 'y';
    end;

end;

if (ansPrint == 'y')
   fprintf(1,'\nPrinting Statistics...');
   print('-dpsc2',namPrint2,fig1) ;   
   fprintf(1,'done');
end;

if (ansPS == 'y')
   fprintf(1,'\nSaving Postscript file...');
   print('-depsc2',fig1,namPS);   
   fprintf(1,'done');
end;

if ishandle(fig1)
    close(fig1);
end;    

clear ansPrint ansPS namPS namPrint namPrint2;
clear bin binTimErr cTol edges hMaxDat ax1 ax2;
clear hiLim hiMaxAx hisMaxTick hiMeanAx hiMeanDat fig1 pos1;
clear hisMeanTick hisTimeErr indMaxErrX indMaxErrY;
clear indMeanErrX indMeanErrY j loLim loMaxAx loMeanAx maxErr maxPlot;
clear n nMaxBin nMeanBin sigp;
%% Calculating the mean flatfield

fprintf(1,'\nCalculating Meanvalue of each pixel...');
fig2 = figure;

sumStack = mean(ImStack,3);
meanFf   = mean2(sumStack);
ssf      = std2(sumStack)/meanFf;
fprintf(1,'done');

ax1 = axes('units','normalized','position',[0.28 .725 0.6 .21]);
   imagesc(sumStack/meanFf);
   axis off;
   set(gca,'clim',[1-2*ssf 1+2*ssf]);
   title('Normalized, mean flatfield image')
   colormap(pilatus);
    
%% Calculating 4 stripes through the flatfield and plot them

uplim1 = round(ydim/2*rand(1)+ydim/2);  
lolim1 = uplim1 - 10;  dimlim1 = uplim1 - lolim1 + 1;
lolim2 = round(ydim/2*rand(1));  
uplim2 = lolim2 + 10;  dimlim2 = uplim2 - lolim2 + 1;

line([0 xdim],[uplim1 uplim1],'Color','r');
line([0 xdim],[uplim2 uplim2],'Color','b');
line([0 xdim],[lolim1 lolim1],'Color','r');
line([0 xdim],[lolim2 lolim2],'Color','b');

lelim1 = round(xdim/2*rand(1));         
rilim1 = lelim1 + 10;  dimlim3 = rilim1 - lelim1 + 1;    
rilim2 = round(xdim/2*rand(1)+xdim/2);
lelim2 = rilim2 - 10;  dimlim4 = rilim2 - lelim2 + 1;
line([lelim1 lelim1],[0 ydim],'Color','c');
line([lelim2 lelim2],[0 ydim],'Color','m');
line([rilim1 rilim1],[0 ydim],'Color','c');
line([rilim2 rilim2],[0 ydim],'Color','m');

axes('units','normalized','position',[0.9 .725 0.04 .21]);
   handle = gca;
   colorbar(handle,'peer',ax1);

axes('units','normalized','position',[0.28 0.575 0.6 0.125]);
   plot(sum(sumStack(lolim1:uplim1,:)/meanFf/(dimlim1),1),'Color','r');
   hold on;
   plot(sum(sumStack(lolim2:uplim2,:)/meanFf/(dimlim2),1),'Color','b');
   xlim([0 xdim]);
   ylim([0.75 1.25]);

axes('units','normalized','position',[0.06 0.725 0.16 0.21]);
   plot(sum(sumStack(:,lelim1:rilim1)/meanFf/(dimlim3),2),'Color','c');
   hold on;
   plot(sum(sumStack(:,lelim2:rilim2)/meanFf/(dimlim4),2),'Color','m');
   xlim([0 ydim]);
   ylim([0.75 1.25]);
   view(90,90);
%% Calculating a corrected flatfield and plot it
clear ax
fprintf(1,'\nCalculating Corrected...');

corrStack = ImStack(:,:,round(zdim*(rand(1)))).*Im;
meanCf    = mean2(corrStack);
fprintf(1,'done.');

ax1 = axes('units','normalized','position',[0.28 0.225 0.6 .21]);
   imagesc(corrStack/meanCf);
   axis off;
   set(gca,'clim',[1-2*ssf 1+2*ssf]);
   title('A flatfield-corrected, randomly chosen flatfield');
   colormap(pilatus);

line([0 xdim],[uplim1 uplim1],'Color','r');
line([0 xdim],[uplim2 uplim2],'Color','b');
line([0 xdim],[lolim1 lolim1],'Color','r');
line([0 xdim],[lolim2 lolim2],'Color','b');

line([lelim1 lelim1],[0 ydim],'Color','c');
line([lelim2 lelim2],[0 ydim],'Color','m');
line([rilim1 rilim1],[0 ydim],'Color','c');
line([rilim2 rilim2],[0 ydim],'Color','m');

axes('units','normalized','position',[0.9 .225 0.04 .21]);
   handle = gca;
   colorbar(handle,'peer',ax1);

axes('units','normalized','position',[0.28 0.075 0.6 0.125]);
   plot(sum(corrStack(lolim1:uplim1,:)/meanCf/(dimlim1),1),'Color','r');
   hold on;
   plot(sum(corrStack(lolim2:uplim2,:)/meanCf/(dimlim1),1),'Color','b');
   xlim([0 xdim]);
   ylim([0.75 1.25]);

axes('units','normalized','position',[0.06 0.225 0.16 0.21]);
   plot(sum(corrStack(:,lelim1:rilim1)/meanCf/(dimlim3),2),'Color','c');
   hold on;
   plot(sum(corrStack(:,lelim2:rilim2)/meanCf/(dimlim4),2),'Color','m');
   xlim([0 ydim]);
   ylim([0.75 1.25]);
   view(90,90);
   
clear ydim xdim zdim
%% Printing and saving dialog for the flatfield

reply = 'n';
ansPrint='';
printerNam = printer;
PSnam      = 'flatfield.eps';

while (reply ~='y')

    ansPrint = input('\nDo you want to print the Figure? [y] :','s');
    if isempty(ansPrint)
        ansPrint = 'y';
    end;

    if (ansPrint=='y')
        namPrint = input(...
            ['What is the name of the Printer? [',printerNam,'] : '],'s');
        if isempty(namPrint)
            namPrint=printerNam;
        end;
        namPrint2=strcat('-P',namPrint);
    end;

    ansPS = input('Do you want to save the Figure as an eps-file? [y] :','s');
    if isempty(ansPS)
        ansPS = 'y';
    end;

    if (ansPS=='y')
        namPS = input( ...
            ['How should the eps-file be called? [',PSnam,'] : '],'s');
        if isempty(namPS)
            namPS=PSnam;
        end;
    end;

    reply = input('Are all inputs OK? y/n [y]: ', 's');
    if isempty(reply)
        reply = 'y';
    end;

end;

if (ansPrint == 'y')
   fprintf(1,'\nPrinting Flatfield...');
   print('-dpsc2',namPrint2,fig2) ;   
   fprintf(1,'done');
end;

if (ansPS == 'y')
   fprintf(1,'\nSaving Postscriptfile...');
   print('-depsc2',fig2,namPS);   
   fprintf(1,'done');
end;

clear ansPrint;
clear ansPS;
clear namPS;
clear namPrint;
clear namPrint2;

if ishandle(fig2)
    close(fig2);
end;

fprintf(1,'\n');

end % of main function

%% subfunction histsummask
%           -----------

% HISTSUMMASK - gets the histmask from the summed histmasks of each and
% every ff image taken. 
%
%
%   [MaskImgOut, <SummedMask>, <nobadpix>] 
%               = histsummask(ImStackIn, thresh, tolerance, <MaskImgIn>)
%
%   MaskImgOut: Maskimg based on tolerance factor operating on SummedMask.
%               Maskimg is 1 for all pixels above "tolerance", 0 otherwise
%   SummedMask: The sum of the n individual masks for each image. Each 
%               image mask contains 1s where pixels that lie within the 
%               limits set by vector "thresh", and zeros otherwise. 
%               Therefore SummedMask is an array that contains integer 
%               values between 0 and n (the number of ff images recorded,
%               and therefore also the size of the 3rd dimension of
%               ImStackIn.
%   nobadpix  : Returns the number of bad pixels in MaskImgOut
%   ImStackIn : Input 3-D stack of images to mask 
%   thresh    : Vector containing both the lower (low) and higher (high)
%               thresholds, thresh = [low high]. Note: pixel intensities
%               identically equal to either low or high are accepted. 
%   tolerance : integer number less than or equal to n, defining the 
%               minimum number of times any given pixel has to lie within
%               the thresholds given by thresh. e.g., for a maximum number
%               of bad pixel events of 10% and 100 images in the stack, set
%               tolerance to 90.
%   MaskImgIn : Optional input mask to select only a part of your 
%               ImStackIn, e.g. to work on a specific region of interest.
%
%   <argument> depicts an optional argument

function [MaskImgOut, SummedMask, nobadpix] = ... 
  histsummask (ImStackIn, thresh, tolerance)

    [ydim xdim zdim]=size(ImStackIn);    zdim;
    MaskImgOut = zeros ([ydim xdim]);
    % sum the zdim histograms of ImStackIn to create SummedMask
    a = histmask (ImStackIn, thresh);
    SummedMask = sum(a,3);

    % search for occasional hot pixels using the tolerance parameter
    indices = (SummedMask >= tolerance);
    MaskImgOut(indices) = 1; % sets bad pixels to zero
    % make Maskimg with ones for all nonzero elements, count nonzero elements
    % in Maskimg
    nobadpix = ydim * xdim - nnz(MaskImgOut);

end % of histsummask

%% subfunction poisson
%           -----------

% POISSON - calculates a theoretical poisson distribution and delivers the 
%           value at position x
%
%   out = poisson(x,lambda)
%   
%   out    : the value at position x of the poissonian
%
%   x      : the position of which you want to habe the value of the
%            poissonian
%   lambda : the mean-value of the distribution

function out = poisson(x,lambda)
   
   if lambda > 100
       out = 1./sqrt(2*pi*x).*(lambda./x).^x.*exp(x-lambda);
   else
       out = lambda.^x ./ gamma(x+1).* exp(-lambda);
   end;

end