function [ppspotdat,ppstdxydat] = ppallspots(immat,spotmat,searchmat,width,imcut)
% Function PPALLSPOTS.M
%
%    [ppspotdat,ppstdxydat] = ppallspots(immat,spotmat,searchmat,width,imcut)
%
%    PPALLSPOTS finds the weighted centres for all the spots within a
%    Pepperpot image IMMAT.  It then calculates the horizontal and vertical
%    offset of the found spots, in [number of pixels], from the expected
%    hole position.  PPALLSPOTS uses the function PPSPOT to find the
%    weighted centre of each light spot, looping over the function 4 times
%    to determine the exact centre of the image spot.  The image matrix
%    IMMAT must be a grayscale image.
%
%    PPALLSPOTS also requires a spot matrix, SPOTMAT, that defines the
%    expected hole positions (in X and Y) of each hole, to allow PPALLSPOTS
%    to calculate the difference.  In addition, a search matrix, SEARCHMAT,
%    can be defined, that starts the search for the spot position from a
%    pixel location other than that defined in SPOTMAT, since using SPOTMAT
%    sometimes results in the wrong spot being found.  The width of the
%    search area, WIDTH, and an image intensity low cut, IMCUT, can also be
%    specified.
%
%    [PPSPOTDAT] = PPALLSPOTS(IMMAT,SPOTMAT) finds the centre of each of
%    the spots in the grayscale image IMMAT, and their offset from the
%    predicted positions, based on the predicted hole positions in SPOTMAT. 
%    The dimensions/arrangement of the matrices are as follows:
%
%        IMMAT - An [A x B] matrix of grayscale image intensity (usually
%                2048 x 2048).
%
%        SPOTMAT - An [M x N x 2] matrix of horizontal and vertical hole
%                  positions, in [number of pixels]: the grid must be M
%                  holes by N holes (usually 25 x 25); [M,N,1] gives the
%                  horizontal pixel location of the hole (M,N), and [M,N,2]
%                  the vertical pixel location.
%
%        PPSPOTDAT - An [M x N x 5] matrix, where (M,N) corresponds to the
%                    hole position as before; [M,N,1:2] gives the positions
%                    of each light spot in absolute number of pixels;
%                    [M,N,3:4] gives the horizontal and vertical pixel
%                    offset for each spot from the predicted hole position
%                    given in SPOTMAT; [M,N,5] is the total light intensity
%                    for that hole.
%
%    [PPSPOTDAT,PPSTDXYDAT] = PPALLSPOTS(IMMAT,SPOTMAT,SEARCHMAT) also
%    outputs PPSTDXYDAT, an [M x N x 2] matrix, with standard deviation
%    information for each spot.  [M,N,1] is the horizontal and [M,N,2] is
%    the vertical.
%
%    [PPSPOTDAT] = PPALLSPOTS(IMMAT,SPOTMAT,SEARCHMAT) specifies the search
%    matrix, SEARCHMAT: this allows the search for the spot to be specified
%    more accurately than just using SPOTMAT.  SEARCHMAT must be the same
%    size as SPOTMAT, since they refer to the same hole arrangement.  The
%    default value is SEARCHMAT = SPOTMAT.
%
%    [PPSPOTDAT] = PPALLSPOTS(IMMAT,SPOTMAT,SEARCHMAT,WIDTH) specifies the
%    width of the search area that is passed to PPSPOT in which to sum up
%    all the image data to find the weighted centre.  The default value is
%    WIDTH = 51.
%
%    [PPSPOTDAT] = PPALLSPOTS(IMMAT,SPOTMAT,SEARCHMAT,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.  The default value is IMCUT = 0.
%

    if nargin < 2
        error('Need at least 2 input variables: [ppspotdat,ppstdxydat] = ppallspots(immat,spotmat,searchmat,width,imcut)') ;
    elseif nargin > 5
        error('Can only specify 5 input variables: [ppspotdat,ppstdxydat] = ppallspots(immat,spotmat,searchmat,width,imcut)') ;
    elseif nargout < 1
        error('Must specify at least 1 output variable: [ppspotdat] = ppallspots(immat,spotmat,searchmat,width,imcut)') ;
    elseif nargout > 2
        error('Can only specify 2 output variables: [ppspotdat,ppstdxydat] = ppallspots(immat,spotmat,searchmat,width,imcut)') ;
    end

    if nargin < 5
        imcut = 0 ;
    end
    if nargin < 4
        width = 51 ;
    end
    if ( nargin < 3 )
        searchmat = spotmat ;
    end
    if isempty(searchmat)
        searchmat = spotmat ;
    end

%    Make sure spotmat and searchmat are the right size

    spmsize = size(spotmat) ;
    semsize = size(searchmat) ;

    if ( length(spmsize) ~= 3 ) || ( spmsize(3) ~= 2 )
        error('spotmat matrix must be an [M x N x 2] matrix for an [M x N] grid of holes') ;
    elseif ( ( sum( spmsize == semsize ) ) ~= 3 )
        error('spotmat and searchmat matrices must be the same size: [M x N x 2]') ;
    end

%    Use PPSPOT to find light spots; loop over each spot 4 times

    ppspotdat = zeros(spmsize(1),spmsize(2),5) ;
    ppstdxydat = zeros(spmsize(1),spmsize(2),2) ;

    for i = 1:spmsize(1)
        for j = 1:spmsize(2)

            spot = [spotmat(i,j,1) spotmat(i,j,2)] ;
            spotb = [searchmat(i,j,1) searchmat(i,j,2)] ;

            [cent,stdxy,spint] = ppspot(immat,spotb,width,imcut) ;
            [cent2,stdxy,spint] = ppspot(immat,cent,width,imcut) ;
            [cent3,stdxy,spint] = ppspot(immat,cent2,width,imcut) ;
            [cent4,stdxy,spint] = ppspot(immat,cent3,width,imcut) ;

            if spint <= 0
                cent4 = spot ;
            end

            ppspotdat(i,j,1:2) = cent4(1:2) ;
            ppspotdat(i,j,3) = ppspotdat(i,j,1) - spot(1) ;
            ppspotdat(i,j,4) = ppspotdat(i,j,2) - spot(2) ;
            ppspotdat(i,j,5) = spint ;
            ppstdxydat(i,j,1:2) = stdxy ;

        end
    end

    return