function [cent,stdxy,summat] = ppcentre(immat,imcut)
% Function PPCENTRE.M
%
%    [cent,stdxy,summat] = ppcentre(immat,imcut)
%
%    PPCENTRE finds the weighted centre of an image, based on the pixel
%    intensities.  The image matrix IMMAT must be a grayscale image.
%
%    PPCENTRE(IMMAT) finds the centre of the image matrix IMMAT and
%    displays the pixel location of that centre.
%
%    PPCENTRE(IMMAT,IMCUT) also applies a cut to the image data, removing
%    pixel intensities below the intensity level IMCUT by subtracting the
%    grayscale value IMCUT from each pixel value, then making the minimum
%    value of any pixel zero.
%
%    [CENT] = PPCENTRE(IMMAT,IMCUT) outputs the central pixel location to
%    the variable CENT.  The central pixel location is not rounded off
%    since CENT is a weighted average.
%
%    [CENT,STDXY] = PPCENTRE(IMMAT,IMCUT) also outputs the horizontal and
%    vertical standard deviations to the variable STDXY.
%
%    [CENT,STDXY,SUMMAT] = PPCENTRE(IMMAT,IMCUT) also outputs the Spot
%    Intensity ie. the sum of all the pixel intensities in the specified
%    image, to the variable SUMMAT.
%

    if nargin < 1
        error('Need at least 1 input variable: [cent,stdxy,summat] = ppcentre(immat,imcut)') ;
    elseif nargin > 2
        error('Can only specify 2 input variables: [cent,stdxy,summat] = ppcentre(immat,imcut)') ;
    elseif nargout > 3
        error('Can only specify 3 output variable: [cent,stdxy,summat] = ppcentre(immat,imcut)') ;
    end

    if nargin == 1
        imcut = 0 ;
    end

%    Apply cut to image intensity

    cutmat = ( immat > imcut ) ;
    immat = ( immat - imcut ).*cutmat ;
    imsize = size(immat) ;

%    Find image centre

    cent = [(imsize(2) - 1)./2 (imsize(1) - 1)./2] ;
    stdxy = [0 0] ;

    summat = sum(sum(immat)) ;

    if summat > 0

        xidx = [1:imsize(2)] ;
        yidx = [1:imsize(1)] ;

        xpix = ( sum( (sum(immat)).*xidx ) )./summat ;
        ypix = ( sum( (sum(immat')).*yidx ) )./summat ;

        cent = [xpix ypix] ;

%    Calculate standard deviations

        xvar = sum(((sum(immat)).*((xidx-xpix).^2)))./summat ;
        yvar = sum(((sum(immat')).*((yidx-ypix).^2)))./summat ;

        stdxy = [sqrt(xvar) sqrt(yvar)] ;
%       summat = sum(sum(cutmat)) ;

    end

    if nargout == 0
        disp(' ')
        disp(['Image centre is at  ' num2str(cent) ' pixels.']) ;
        disp(' ')
    end

    return