function [ppim,stdm,pphist,ppminmax] = pepperplotfun(imfile,ims,plt,bitdiv,beamview)
% Function PEPPERPLOTFUN.M
%
%    [ppim,stdm,pphist,ppminmax] = pepperplotfun(imfile,ims,plt,bitdiv,beamview)
%
%    PEPPERPLOT loads image data from the pepperpot emittance scanner and
%    plots histograms of the image data.  The data is loaded from multiple
%    image TIFF files.
%
%    [PPIM] = PEPPERPLOTFUN(IMFILE) loads the first image in the image file
%    IMFILE, displays the image, rescales the image to enhance the features
%    and displays it and calculates and displays a histogram of the pixel
%    intensity information.
%
%    [PPIM] = PEPPERPLOTFUN(IMFILE,N) loads the Nth image from the image
%    file.  If the number of images is less than N, the last image in the
%    file is loaded.  N can be a vector of image numbers ie. [4 6 10] loads
%    the 4th, 6th and 10th image.
%
%    [PPIM] = PEPPERPLOTFUN(IMFILE,0) loads all the images from the image
%    file and calculates an average and a standard deviation.
%
%    [PPIM] = PEPPERPLOTFUN(IMFILE,N,[PLT]) allows control of the image and
%    histogram plotting.  PLT needs to be a 2-element variable: PLT = [A B].
%    If A == 1, the main image is plotted; if A == 2, the main image is
%    plotted with scaled intensities.  If B == 1, the histogram of
%    grayscale values and the standard deviation image are plotted; if B == 2,
%    the histogram of grayscale values in plotted on its own plot.  Both
%    variables must be specified ie. length(PLT) = 2.
%
%    [PPIM] = PEPPERPLOTFUN(IMFILE,N,[PLT],BITDIV) allows control over the
%    bit depth of the tiff image.  Normal tiff images are 16-bit, but those
%    output by the PCO camera are 14-bit and need to be divided by 4 to get
%    the correct bit values.  For 16-bit images, BITDIV = 1; for 14-bit PCO
%    images, BITDIV = 4: this is the default value.
%
%    [PPIM] = PEPPERPLOTFUN(IMFILE,N,[PLT],BITDIV,BEAMVIEW) selects the
%    orientation of the image with the BEAMVIEW variable.  Emittance data
%    should be calculated using beam's-eye-view data; camera's-eye-view
%    data is easier for viewing.  Setting BEAMVIEW = ('CAM') loads and
%    displays the image as it was viewed by the camera; BEAMVIEW = ('BEAM')
%    flips the image left-to-right, using FLIPLR(PPIM), to show the
%    beam's-eye-view image.  BEAMVIEW = ('CAM') is the default value.
%
%    [PPIM,STDM] = PEPPERPLOTFUN(IMFILE,0) outputs the standard deviations
%    for the image and displays this data as a scaled image.
%
%    [PPIM,STDM,PPHIST] = PEPPERPLOTFUN(IMFILE,0) outputs histogram data
%    for each image, using one bin for each grayscale level.
%
%    [PPIM,STDM,PPHIST,PPMINMAX] = PEPPERPLOTFUN(IMFILE,0) outputs the
%    minimum and maximum grayscale level for the image: 
%
%            PPMINAX = [MINGS MAXGS]
%
%    The function requires the following input variables:
%
%    IMFILE - The data/image file with the image information, including
%              the file extension and the path (absolute or relative) -
%              must be a string.
%    IMS - The image number(s) within the image file to be read in
%          (multi-image file); the default is 0.  Set this variable to 0 to
%          read out all images.
%    PLT - Variable to specify whether data is to be plotted or not.
%          Default is [2 1] (all data is plotted with main image scaled).
%

    if nargin < 1
        error('Need at least 1 input variable: [ppim,stdm,pphist,ppminmax] = pepperplotfun(imfile,ims,plt)') ;
    elseif nargin > 5
        error('Too many input variables: [ppim,stdm,pphist,ppminmax] = pepperplotfun(imfile,ims,plt,bitdiv,beamview)') ;
    elseif nargout > 4
        error('You can only specify 4 output variables: [ppim,stdm,pphist,ppminmax] = pepperplotfun(imfile,ims,plt)') ;
    elseif ( nargin == 3 ) & ( length(plt) ~= 2 )
        error('You must specify both plotting variables: pepperplotfun(imfile,ims,[2 1])') ;
    elseif ~ischar(imfile)
        error('Image filename must be a string eg. [''imfile.tif'']') ;
    end

    if nargin < 5 || ( isempty(beamview) )
        beamview = ('cam') ;
    end
    if ( nargin < 4 ) || ( isempty(bitdiv) )
        bitdiv = 4 ;
    end
    if nargin < 3 || ( isempty(plt) )
        plt = [2 1] ;
    end
    if nargin < 2
        ims = 0 ;
    end

    imdata = imfinfo(imfile,'tif') ;
    width = imdata(1).Width ;
    height = imdata(1).Height ;
    imnum = length(imdata) ;
    ims = fix(ims) ;

    if ( ims(end) > imnum ) | ( length(ims) > imnum )
        error('Too many images or too large an image number specified in IMS variable') ;
    end

    if ims > imnum
        ims = imnum ;
    elseif ims <= 0
        ims = [1:1:imnum] ;
    end

    imloop = length(ims) ;    

%    Create necessary variables

    ppimm = zeros(height,width) ; ppimv = zeros(height,width) ;
    binmids = [1:1:16384] ;

%    Load image file

    if ( length(ims) > 1 )

        disp(' ') ;
        disp(['Loading ' num2str(imloop) ' images from image file ' imfile ' ... ']) ;
        disp(' ') ;

        for i = 1:imloop
            disp(['Loading image ' num2str(ims(i)) ' ... ']) ;
            [ppimu(:,:)] = imread(imfile,'tif',ims(i)) ;
            ppimd = double(ppimu)./bitdiv ;
            ppimm = ppimm + ppimd ;
            ppimv = ppimv + ppimd.^2 ;
        end

        ppim = ppimm./imloop ;
        stdm = sqrt( ( ppimv./imloop ) - ppim.^2 ) ;

        disp(' ') ;

    elseif ( ( ims > 0 ) | ( imnum == 1 ) )
        
        disp(' ') ;
        disp(['Loading image ' num2str(ims) ' from image file ' imfile ' ... ']) ;
        disp(' ') ;

        [ppimu(:,:)] = imread(imfile,'tif',ims) ;

        ppim = double(ppimu)./bitdiv ;
        stdm = ppimv ;

    end

    if strcmp(beamview,'beam')
        ppim = fliplr(ppim) ;
        stdm = fliplr(stdm) ;
    end

%    imsize = size(ppim) ;
    ppimres = reshape(ppim,[],1) ;
    maxgs = squeeze(max(max(ppim))) ;
    mings = squeeze(min(min(ppim))) ;
    maxgstd = squeeze(max(max(stdm))) ;
    mingstd = squeeze(min(min(stdm))) ;
    ppplot = squeeze(ppimres) ;

    ppminmax = [mings maxgs] ;

    [pphist] = hist(ppplot,binmids) ;

%    Create new figure, rescale and plot image

    if plt(1) == 1

        figure, imagesc(ppim,[1 16384]), colormap(gray), axis image ;
        title(['\bf Raw Image Data']) ;

    elseif plt(1) == 2

        figure, imagesc(ppim,[mings maxgs]), colormap(gray), axis image ;
        title(['\bf Scaled Image Data']) ;

    end

    if plt(2) == 1

        figh = figure ;
        screen = get(0,'Screensize') ;
        figpos = screen ; figpos(4) = figpos(4).*0.9 ;
%        set(figh,'Position',figpos) ;

        figure(figh), subplot(2,3,1), imagesc(ppim(:,:),[1 16384]), colormap(gray), axis image ;
        title(['\bf Raw Image Data']) ;
        figure(figh), subplot(2,3,4), imagesc(ppim(:,:),[mings maxgs]), colormap(gray), axis image ;
        title(['\bf Scaled Image Data']) ;

%    Create histogram of image intensity data and plot

        figure(figh), subplot(2,3,[2 6]), bar(binmids,pphist,'k') ;
        histax = axis ;
        histax(1) = min(mings) ; histax(2) = max(maxgs) ;
        if histax(1) >= histax(2)
            histax(1) = histax(1) - 25 ;
            histax(2) = histax(1) + 50 ;
        end
        axis(histax) ;
        xlabel(['16-bit light intensity per pixel (counts)']), ylabel(['Number of pixels']), ...
            title(['\bf Light intensity distribution for pepperpot image']) ;

%    Display standard deviations as image

        if ( length(ims) > 1 )

            fighh = figure ;
            imagesc(stdm(:,:),[mingstd maxgstd]), colormap(gray), axis image ;
            title(['\bf Standard deviations']) ;

        end

    elseif plt(2) == 2

%    Create histogram of image intensity data and plot

        figh = figure ;
        bar(binmids,pphist,'k') ;
        histax = axis ;
        histax(1) = min(mings) ; histax(2) = max(maxgs) ;
        if histax(1) >= histax(2)
            histax(1) = histax(1) - 25 ;
            histax(2) = histax(1) + 50 ;
        end
        axis(histax) ;
        xlabel(['14-bit light intensity per pixel (counts)']), ylabel(['Number of pixels']), ...
            title(['\bf Light intensity distribution for pepperpot image']) ;
    end

    return




