function [spotim,gaussim] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)
%
% function [spotim,gaussim] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)
%
%    GAUSSPIXSPOT - generate a black and white pixel image of a Gaussian
%    light spot.
%
%    [spotim] = gausspixspot(imsize,sigmapix)
%
%    [spotim,gaussim] = gausspixspot(imsize,sigmapix)
%
%    [...] = gausspixspot(imsize,sigmapix,sigmacut)
%    [...] = gausspixspot(imsize,sigmapix,sigmacut,maxint)
%    [...] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset)
%    [...] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)
%
%    GAUSSPIXSPOT generates a pixel image of a light spot that obeys a
%    Gaussian distribution in both directions.  This simulates the light
%    that falls on a pixel, or array of pixels, in a CCD camera image from
%    a point source: the light that spreads out from this source normally
%    obeys a Gaussian distribution.
%
%    To simulate the image of such a spot, it is necessary to integrate a
%    Gaussian in 2 directions - both X and Y - in order to calculate the
%    amount of light that is collected by each pixel.  This requires the
%    calculation of 2 error functions, ERF(X) and ERF(Y).
%
%    [spotim] = gausspixspot(imsize,sigmapix) - output a Gaussian spot
%    image to the image array SPOTIM.  The size of the image, in number of
%    pixels, is given by IMSIZE: this can either be a single value or a 
%    [1 x 2] array.  If a single value is specified, the image is assumed
%    to be symmetrical; if a [1 x 2] array is used, this specifies the 
%    [X Y] size of the image.  The standard deviation of the Gaussian, in
%    number of pixels, is given by SIGMAPIX: as with IMSIZE, [SIGMA] gives
%    both the X & Y sizes of the distribution, while [SIGMAX SIGMAY] gives
%    the X and Y sizes independently.
%
%    [spotim,gaussim] = gausspixspot(imsize,sigmapix) - also outputs a
%    Gaussian spot image to the image array GAUSSIM.  Instead of
%    integrating over each pixel area, GAUSSIM is produced straight from
%    the magnitude of the initial Gaussian distribution at the centre of
%    each pixel.  For large numbers of pixels, SPOTIM and GAUSSIM are a
%    good approximation.  GAUSSIM is mainly used for comparative purposes.
%
%    [...] = gausspixspot(imsize,sigmapix,sigmacut) - also specifies the
%    number of sigma at which to cut off the distribution with SIGMACUT.
%    SIGMACUT is only approximate, since it isn't possible to calculate
%    partially filled pixels: any pixels that partially fall inside
%    SIGMACUT will be calculated.  As before, specifying [SIGMACUT] assumes
%    the same cutoff in both X and Y, while [SIGMACUTX SIGMACUTY] specifies
%    the cutoff values for X and Y independently.  Setting SIGMACUT = 0
%    means that no cutoff is applied: this is the default value.  
%
%    Note that, because of the fact that the distributions may be different
%    in X and Y, SIGMACUT will result in a square distribution, not a
%    cutoff at a given radius.
%
%    [...] = gausspixspot(imsize,sigmapix,sigmacut,maxint) - set the
%    maximum intensity of the peak of the distribution to be MAXINT.  For 
%    example, for a fully filled 16-bit grayscale image, MAXINT = 65,536.
%    The default valye is MAXINT = 1: pixel intensities all fall in the
%    range [0:1].
%
%    [...] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset) - specify
%    the offset, in number of pixels, of the centre of the distribution
%    from the centre of the image.  While this can be used to shift the
%    distribution a large distance from the centre of the image, the main
%    use is to move the distribution by fractions of a pixel so that it is
%    no longer symmetric in its pixel distribution.  The default value is
%    OFFSET = [0,0]: image is centred.
%
%    Note that vertical offset values should be the negative of what you
%    may expect, since pixel values start from the top-left, not the
%    bottom-left, so the vertical axis runs from the top down.
%
%    [spotim] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)
%    - round off the pixel intensities to integer values by specifying
%    ROUNDUP = 'round'.  This approximates real images, where pixels can
%    only have integer grayscale values: this is only really useful when
%    MAXINT is some large number.
%
%    BACKGROUND
%
%    In 1 dimension, a Gaussian distribution is given by:
%
%        y = A e^-[ (x-mu)^2 / 2*(sigma^2) ]
%
%    ...where mu is the mean, sigma is the standard deviation and A is an
%    amplitude.  In a fully normalised Gaussian distribution, 
%    A = 1/(sqrt(2*pi)).  The integral of the Gaussian is given by the
%    Error Function, ERF(X), where:
%
%        1/2*erf( (x-mu)/(sigma*sqrt(2)) )
%
%    ...is the integral of the Gaussian distribution between 0 and x.
%
%    In 2 dimensions, the Gaussian distribution is given by:
%
%        z = A e^-[(x-mux)^2 / 2*(sigmax^2) + (y-muy)^2 / 2*(sigmay^2)]
%          = A e^-[(x-mux)^2 / 2*(sigmax^2)] e^-[(y-muy)^2 / 2*(sigmay^2)]
%
%    ...where mux/muy are the means in the X/Y directions, sigmax/sigmay is
%    the standard deviation in the X/Y directions, and A is the amplitude
%    as before.  As such, the integral of this distribution is given by:
%
%        1/4*[erf((x-mux)/(sigmax*sqrt(2)))]*[erf((y-muy)/(sigmay*sqrt(2)))]
%
%    As before, X and Y must be between 0 and infinity.  As such, the
%    calculation is made in 4 parts: for positive and negative X and Y.

% File is the property of Simon Jolly, Imperial College London.  File
% released under the GNU public license.
%
% File history
%
%   30-Sep-2010 S. Jolly
%       Creation of fake pepperpot calibration image for new 4 light spot
%       calibration system.
%
%======================================================================


    if nargin < 2
        error('Need at least 2 input variables: [spotim,gaussim] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)') ;
    elseif nargin > 6
        error('Can only specify 6 input variables: [spotim,gaussim] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)') ;
    elseif nargout > 2
        error('Can only specify 2 output variables: [spotim,gaussim] = gausspixspot(imsize,sigmapix,sigmacut,maxint,offset,roundup)') ;
    end

    if nargin < 6 || isempty(roundup)
        roundup = 0 ;
    elseif strncmpi(roundup,'round',5)
        roundup = 1 ;
    else
        roundup = 0 ;
    end

    if nargin < 5 || isempty(offset)
        offset = [0 0] ;
    elseif numel(offset) == 1
        offset = [offset -offset] ;
    elseif numel(offset) > 2
        error('Input variable OFFSET can only contain 1 or 2 values') ;
    end

    if nargin < 4 || isempty(maxint)
        maxint = 1 ;
    elseif numel(maxint) > 1
        error('Input variable MAXINT must contain a single value') ;
    end

    if nargin < 3 || isempty(sigmacut)
        sigmacut = [0 0] ;
    elseif numel(sigmacut) == 1
        sigmacut = [sigmacut sigmacut] ;
    elseif numel(sigmacut) > 2
        error('Input variable SIGMACUT can only contain 1 or 2 values') ;
    end

    if numel(sigmapix) == 1
        sigmapix = [sigmapix sigmapix] ;
    elseif numel(sigmacut) > 2
        error('Input variable SIGMAPIX can only contain 1 or 2 values') ;
    end

    if numel(imsize) == 1
        imsize = [imsize imsize] ;
    elseif numel(imsize) > 2
        error('Input variable IMSIZE can only contain 1 or 2 values') ;
    end

%    Create variables for image

    initim = zeros(imsize) ;
    spotim = initim ;
    gaussim = initim ;
    intim = initim ;
    xpixedges = [0:1:imsize(1)] ;
    ypixedges = [0:1:imsize(2)] ;
    sigmax = sigmapix(1) ;
    sigmay = sigmapix(2) ;

    midpoint = imsize./2 ;
    midoff = midpoint + offset ;
    mux = midoff(1) ;
    muy = midoff(2) ;

%    Calculate intensities by cheating

    xpixmidvals = xpixedges(2:end) - 0.5 ;
    ypixmidvals = ypixedges(2:end) - 0.5 ;
    [xpixmids, ypixmids] = meshgrid(xpixmidvals, ypixmidvals) ;

    gaussim = ( exp( -( (xpixmids - mux).^2 )./(2.*sigmax.^2) ) ).*( exp( -( (ypixmids - muy).^2 )./(2.*sigmay.^2) ) ) ;
    maxgim = max(gaussim(:)) ;
    gaussim = gaussim.*(maxint./maxgim) ;

%    Calculate intensities by integration

    xstartvals = xpixedges(1:end-1) ;
    ystartvals = ypixedges(1:end-1) ;
    xendvals = xpixedges(2:end) ;
    yendvals = ypixedges(2:end) ;
    [xstarts, ystarts] = meshgrid(xstartvals, ystartvals) ;
    [xends, yends] = meshgrid(xendvals, yendvals) ;

    intimfull = 0.25.*erf( (xends - mux)./(sigmax.*sqrt(2)) ).*erf( (yends - muy)./(sigmay.*sqrt(2)) ) ;
    intimshort = 0.25.*erf( (xends - mux)./(sigmax.*sqrt(2)) ).*erf( (ystarts - muy)./(sigmay.*sqrt(2)) ) ;
    intimprefull = 0.25.*erf( (xstarts - mux)./(sigmax.*sqrt(2)) ).*erf( (yends - muy)./(sigmay.*sqrt(2)) ) ;
    intimpreshort = 0.25.*erf( (xstarts - mux)./(sigmax.*sqrt(2)) ).*erf( (ystarts - muy)./(sigmay.*sqrt(2)) ) ;
    intimslice = intimfull - intimshort ;
    intimpre = intimprefull - intimpreshort ;
    spotim = intimslice - intimpre ;
    maxim = max(spotim(:)) ;
    spotim = spotim.*(maxint./maxim) ;

    plusonesx = initim ;
    plusonesy = initim ;

    if sigmacut(1) > 0

        lowx = mux - (sigmacut(1).*sigmax) ;
        highx = mux + (sigmacut(1).*sigmax) ;
        indsx =  xstarts >= lowx & xends <= highx  ;
        plusonesx(indsx) = 1 ;

    else

        plusonesx = 1 ;

    end

    if sigmacut(2) > 0

        lowy = muy - (sigmacut(2).*sigmay) ;
        highy = muy + (sigmacut(2).*sigmay) ;
        indsy =  ystarts >= lowy & yends <= highy  ;
        plusonesy(indsy) = 1 ;

    else

        plusonesy = 1 ;

    end

    spotim = spotim.*plusonesx.*plusonesy ;
    gaussim = gaussim.*plusonesx.*plusonesy ;

    if roundup
        spotim = round(spotim) ;
        gaussim = round(gaussim(:,:)) ;
    end

    return








