function [cent,stdxy,spint] = ppspot(immat,spot,width,imcut)
% Function PPSPOT.M
%
%    [cent,stdxy,spint] = ppspot(immat,spot,width,imcut)
%
%    PPSPOT finds the weighted centre of a set of pixels within an image.
%    The image matrix IMMAT must be a grayscale image.  PPSPOT calls the
%    function PPCENTRE to calculate the weighted centre of a group of
%    sub-pixels.  This sub-group is selected from the main image matrix by
%    specifying a centre point and a width.  In addition, a cut can be
%    applied to the intensity.
%
%    PPSPOT(IMMAT,SPOT,WIDTH) finds the centre of an area within the image
%    matrix IMMAT.  The centre of the area is defined by SPOT, where 
%    SPOT = [XPIX YPIX] and the width of the area to search (both
%    horizontal and vertical) is defined by WIDTH pixels.  The area is then
%    displayed on screen.  WIDTH can either be a pair of values
%    corresponding to the width and height of the selected region ie.
%    WIDTH = [XSIZE YSIZE], or a single value that defines both the
%    horizontal and vertical size of the search area.
%
%    Since PPSPOT only deals in whole pixels, WIDTH should be an odd
%    number; otherwise the width is rounded up to the nearest odd number.
%    The pixel location SPOT should be a whole number: PPSPOT automatically
%    rounds it to the nearest whole pixel number before creating the
%    sub-image.
%
%    PPSPOT(IMMAT,SPOT,WIDTH,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] = PPSPOT(IMMAT,SPOT,WIDTH,IMCUT) outputs the central pixel
%    location to the variable CENT.  CENT is a weighted average and is not
%    rounded off to the nearest pixel value.
%
%    [CENT,STDXY] = PPSPOT(IMMAT,SPOT,WIDTH,IMCUT) also outputs the
%    horizontal and vertical standard deviations to the variable STDXY.
%
%    [CENT,STDXY,SPINT] = PPSPOT(IMMAT,SPOT,WIDTH,IMCUT) also outputs the
%    Spot Intensity ie. the sum of all the pixel intensities in the
%    specified image, to the variable SPINT.
%

    if nargin < 3
        error('Need at least 3 input variables: [cent,stdxy,spint] = ppspot(immat,spot,width)') ;
    elseif nargin > 4
        error('Can only specify 4 input variables: [cent,stdxy,spint] = ppspot(immat,spot,width,imcut)') ;
    elseif nargout > 3
        error('Can only specify 2 output variables: [cent,stdxy,spint] = ppspot(immat,spot,width,imcut)') ;
    end

    if nargin < 4
        imcut = 0 ;
    end

    imsize = size(immat) ;

%    Select sub-pixel image from main matrix using spot position

    if numel(width) == 2

        xrad = round( ( width(1) - 1)./2 ) ;
        yrad = round( ( width(2) - 1)./2 ) ;

    elseif numel(width) == 1
        
        xrad = round( ( width - 1)./2 ) ;
        yrad = round( ( width - 1)./2 ) ;

    else

        error('WIDTH variable must contain either a single value or a pair of X and Y values.') ;

    end

    spotro = round( spot ) ;
%    spotrad = round( ( width - 1)./2 ) ;
%    disp(['X: ' num2str(spotro(1) - xrad) ' to ' num2str(spotro(1) + xrad) ' pixels']) ;
%    disp(['Y: ' num2str(spotro(2) - yrad) ' to ' num2str(spotro(2) + yrad) ' pixels']) ;

    if spotro(2) <= 0 || spotro(1) <= 0 || spotro(2) > imsize(1) || spotro(1) > imsize(2)

        warning('Spot centre falls outside image') ;

        cent = spot ;
        stdxy = 0 ;
        spint = 0 ;

%    elseif ( ( spotro(2) - yrad ) <= 0 ) | ( ( spotro(1) - xrad ) <= 0 ) | ...
%            ( ( spotro(2) + yrad ) > imsize(1) ) | ( ( spotro(1) + xrad ) > imsize(2) )
%
%        warning('Spot radius falls outside image: using default') ;
%
%        cent = spot ;
%        stdxy = 0 ;
%        spint = 0 ;

    else

        if ( spotro(2) - yrad ) <= 0
            ymin = 1 ;
        else
            ymin = spotro(2) - yrad ;
        end
        if ( spotro(2) + yrad ) > imsize(1)
            ymax = imsize(1) ;
        else
            ymax = spotro(2) + yrad ;
        end
        if ( spotro(1) - xrad ) <= 0
            xmin = 1 ;
        else
            xmin = spotro(1) - xrad ;
        end
        if ( spotro(1) + xrad ) > imsize(2)
            xmax = imsize(2) ;
        else
            xmax = spotro(1) + xrad ;
        end

%        subim = immat(((spotro(2) - spotrad):(spotro(2) + spotrad)), ((spotro(1) - spotrad):(spotro(1) + spotrad))) ;
%        subim = immat(((spotro(2) - yrad):(spotro(2) + yrad)), ((spotro(1) - xrad):(spotro(1) + xrad))) ;
        subim = immat((ymin:ymax), (xmin:xmax)) ;

%    Use PPCENTRE to find image centre and standard deviations

        [spotcent,stdxy,spint] = ppcentre(subim,imcut) ;

        if spint <= 0
            cent = spot ;
        else
%            cent = [(spotcent(1) + spotro(1) - spotrad - 1) (spotcent(2) + spotro(2) - spotrad - 1)] ;
%            cent = [(spotcent(1) + spotro(1) - xrad - 1) (spotcent(2) + spotro(2) - yrad - 1)] ;
            cent = [(spotcent(1) + xmin - 1) (spotcent(2) + ymin - 1)] ;
        end

    end

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

    return