function [midpos,sp2sp,imtilt,tout,bout,lout,rout] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs,pixframe)

% function [midpos,sp2sp,imtilt,tout,bout,lout,rout] = 
%    ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs,pixframe)
%
%    PPSPOTCAL - calculate calibration data for pepperpot from 4-light-spot
%    calibration image.
%
%    ppspotcal(immat)
%
%    [midpos] = ppspotcal(immat)
%    [midpos,sp2sp] = ppspotcal(immat)
%    [midpos,sp2sp,imtilt] = ppspotcal(immat)
%
%    [...] = ppspotcal(immat,spotdist)
%    [...] = ppspotcal(immat,spotdist,cutlev)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs)
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs,pixframe)
%
%    [midpos,sp2sp,imtilt,tout,bout,lout,rout] = ppspotcal(...)
%
%    PPSPOTCAL analyses a pepperpot calibration image and calculates the
%    size, centre and orientation of the pepperpot head.  The calibration
%    image must be a TIFF image, taken of the pepperpot head using the 
%    4-light calibration system.  It is always assumed that the calibration
%    spots sit at the top, bottom, left and right of the image: the
%    spot-finding code WILL NOT WORK if the spots are positioned in the
%    corners.
%
%    ppspotcal(immat) - load the image matrix or the image file IMMAT, find
%    the positions of the 4 calibration light spots and display the
%    resulting calibration image.  IMMAT can either be an [M x N] image
%    array in the Matlab workspace that contains image intensity data, or
%    the location of the calibration image file in the form of a string: if
%    IMMAT is a string, PPSPOTCAL calls PEPPERPLOTFUN to load the image
%    data.  It is assumed that the calibration markings are seperated by a
%    distance of 140 mm top-to-bottom and left-to-right (see SPOTDIST
%    below).
%
%    [midpos] = ppspotcal(immat) - outputs the centre position of the
%    calibration image to the variable MIDPOS.
%
%    [midpos,sp2sp] = ppspotcal(immat) - also outputs the spot-to-spot
%    spacing of the mesh, in pixels, to the variable SP2SP.  This assumes
%    a default hole spacing in the pepperpot intercepting mesh of 3 mm (see
%    HOLEPITCH below).
%
%    [midpos,sp2sp,imtilt] = ppspotcal(immat) - also outputs the tilt of
%    the image in radians to the variable IMTILT.
%
%    [...] = ppspotcal(immat,spotdist) - specify the distance, in MM,
%    between the top-to-bottom and left-to-right calibration spots.  The
%    default is 140 mm.
%
%    [...] = ppspotcal(immat,spotdist,cutlev) - apply a cut to the image,
%    removing CUTLEV counts from every pixel.  Any pixel below 0 is set to
%    0 so that these pixels are ignored by the spot-finding code.  The
%    default value is CUTLEV = 0: no cut applied.
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize) - specify the size,
%    SPOTSIZE, in number of pixels, of the region surrounding the brightest
%    pixel in each region.  The default value is SPOTSIZE = 45.
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut) - removes
%    spots from the calibration calculation that have RMS X/Y distribution
%    that is greater than SPOTRMSCUT pixels.  This helps to remove light
%    spots that will give poor calibration measurements because of extra
%    satellite spots: if the Gaussian peak formed by the light from the
%    calibration spot is too big, or too mishapen, it is assumed that too
%    much light comes from outside the central spot and the spot is
%    discarded.  If SPOTRMSCUT <= 0, no cut is applied: this is the default
%    value.
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch)
%    - specify the pitch of the holes in the intercepting pepperpot mesh
%    with HOLEPITCH, given in MM.  This allows scaling of the SP2SP output
%    variable.  The default value is HOLEPITCH = 3 MM.
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh)
%    - select whether the calibration image should be plotted with CALFIGH.
%    It is possible to use an existing figure number on which to plot the
%    image or specify whether or not the image should be plotted at all.
%    CALFIGH can have the following values:
%
%        CALFIGH = [] (empty matrix) or not specifed: don't plot image
%        CALFIGH = 0: plot image in new figure
%        CALFIGH > 0: assumes that CALFIGH is an existing figure and plots
%            the image in figure(calfigh)
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv)
%    - specify the bit depth of the TIFF image file to be read out with
%    BITDIV: this is only for BIG ENDIAN image files ie. for a 14-bit image
%    that is stored as a 16-bit image, BITDIV = 4.  For all little endian
%    files or true 16-bit images, BITDIV = 1: this is the default value.
%    BITDIV is only used when IMMAT specifies an image file location and is
%    passed directly to PEPPERPLOTFUN that loads the image data.
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs)
%    - set the minimum and maximum grayscale levels of the calibration
%    image for plotting with MINMAXGS: this is only necessary if CALFIGH is
%    set to plot the image.  MINMAXGS must be a 2-element array, specifying
%    a minimum and maximum grayscale value for the plotted calibration
%    image eg. MINMAXGS = [200 5000].  If MINMAXGS is empty or not
%    specified, the default grayscale values are taken from the minimum and
%    maximum levels of the image ie. min(MINMAXGS) = min(IMMAT(:)) and
%    max(MINMAXGS) = max(IMMAT(:)).
%
%    [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs,pixframe)
%    - put a mask around the image PIXFRAME pixels wide.
%
%    [midpos,sp2sp,imtilt,tout,bout,lout,rout] = ppspotcal(...) - output
%    the pixel locations of the top spot, TOUT, bottom spot, BOUT, left
%    spot, LOUT, and the right spot, ROUT.

% File is the property of Simon Jolly, Imperial College London.  File
% released under the GNU public license.
%
% File history
%
%   13-Oct-2010 S. Jolly
%       Analysis of 4-light-spot calibration image and output of image
%       data.
%
%======================================================================


    if nargin < 1
        error('Need at least 1 input variable: [...] = ppspotcal(immat)') ;
    elseif nargin > 10
        error('Can only specify 10 input variables: [...] = ppspotcal(immat,spotdist,cutlev,spotsize,spotrmscut,holepitch,calfigh,bitdiv,minmaxgs,pixframe)') ;
    elseif nargout > 7
        error('Can only specify 7 output variable: [midpos,sp2sp,imtilt,tout,bout,lout,rout] = ppspotcal(immat)') ;
    end

    if nargin < 10 || isempty(pixframe)
        pixmask = 0 ;
        pixframe = 0 ;
    elseif sum(pixframe) > 0
        pixmask = 1 ;
        if length(pixframe) ~= 1 && length(pixframe) ~= 4
            error('PIXRAME variable must be either single value or [1 x 4] array') ;
        end
    else
        pixmask = 0 ;
    end

    if nargin < 8 || isempty(bitdiv)
        bitdiv = 1 ;
    end

    if nargin < 7 || isempty(calfigh)
        doplot = 0 ;
        makeplot = 0 ;
    elseif calfigh == 0
        doplot = 1 ;
        makeplot = 1 ;
    else
        doplot = 1 ;
        makeplot = 0 ;
    end

    if nargin < 6 || isempty(holepitch)
        holepitch = 3 ;
    end

    if nargin < 5 || isempty(spotrmscut)
        spotrmscut = 0 ;
    end

    if nargin < 4 || isempty(spotsize)
        spotsize = 45 ;
    end

    if nargin < 3 || isempty(cutlev)
        cutlev = 0 ;
    end

    if nargin < 2 || isempty(spotdist)
        spotdist = 140 ;
    end

%    Set up variables

    calim = [] ;
    tgood = 0 ; bgood = 0 ; lgood = 0 ; rgood = 0 ;

%    Load pepperpot calibration image

    if ischar(immat) && ( exist(immat,'file') == 2 )
        [calim] = pepperplotfun(immat,0,[0 0],bitdiv,'beam') ;
    elseif isnumeric(immat)
        calim = immat ;
    else
        error('IMMAT is neither a valid filename nor a valid Matlab workspace variable.') ;
    end

%    Apply pixel mask and calculate grayscale range

    if pixmask

        if length(pixframe) == 1
            pixframe = [pixframe pixframe pixframe pixframe] ;
        end

        if pixframe(1) > 0
            calim(1:pixframe(1),:) = 0 ;
        end
        if pixframe(2) > 0
            calim(end-pixframe(2)+1:end,:) = 0 ;
        end
        if pixframe(3) > 0
            calim(:,1:pixframe(3)) = 0 ;
        end
        if pixframe(4) > 0
            calim(:,end-pixframe(4)+1:end) = 0 ;
        end

    end

    if nargin < 9 || isempty(minmaxgs)
        calmin = min(calim(:)) ; 
        calmax = max(calim(:)) ;
    else
        calmin = min(minmaxgs(:)) ; 
        calmax = max(minmaxgs(:)) ;
    end        

%    Cut image into 9 segments and find calibration spots

    imheight = size(calim,1) ;
    imwidth = size(calim,2) ;
    segwidth = round(imwidth./6).*2 - 1 ;
    segheight = round(imheight./6).*2 - 1 ;

%    midpos = [imwidth./2 imheight./2] ;
    midpos = [0 0] ;
    midpred = [0 0] ;

%    xsegs = [1 round(imwidth./3) round(2.*imwidth./3) imwidth] ;
%    ysegs = [1 round(imheight./3) round(2.*imheight./3) imheight] ;

    topcent = [round(imwidth./2) round(imheight./6)] ;
    topspot = topcent ;
    botcent = [round(imwidth./2) round(5.*imheight./6)] ;
    botspot = botcent ;
    leftcent = [round(imwidth./6) round(imheight./2)] ;
    leftspot = leftcent ;
    rightcent = [round(5.*imwidth./6) round(imheight./2)] ;
    rightspot = rightcent ;

%    Find top spot

    [topcent1,stdxy,spint] = ppspot(calim,topcent,[segwidth segheight],cutlev) ;
    topwidth = 2.*topcent1(2) ;
    [topcent2,stdxy,spint] = ppspot(calim,topcent1,topwidth,cutlev) ;
    if spotsize > topwidth
        warning('Using smaller area than specified to find top spot: area falls outside image') ;
        topspotsize = topwidth ;
    else
        topspotsize = spotsize ;
    end
    [topcent3,stdxy,spint] = ppspot(calim,topcent2,topspotsize,cutlev) ;
    [topcent4,stdxyt,spintt] = ppspot(calim,topcent3,topspotsize,cutlev) ;

    if spotrmscut <= 0
        stdxytlim = max(stdxyt) + 1 ;
    else
        stdxytlim = spotrmscut ;
    end

    if spintt > 0 && max(stdxyt) < stdxytlim
        topspot = topcent4 ;
        tgood = 1 ;
    end

%    Find bottom spot

    [botcent1,stdxy,spint] = ppspot(calim,botcent,[segwidth segheight],cutlev) ;
    botwidth = 2.*( imheight - botcent1(2) ) ;
    [botcent2,stdxy,spint] = ppspot(calim,botcent1,botwidth,cutlev) ;
    if spotsize > botwidth
        warning('Using smaller area than specified to find bottom spot: area falls outside image') ;
        botspotsize = botwidth ;
    else
        botspotsize = spotsize ;
    end
    [botcent3,stdxy,spint] = ppspot(calim,botcent2,botspotsize,cutlev) ;
    [botcent4,stdxyb,spintb] = ppspot(calim,botcent3,botspotsize,cutlev) ;

    if spotrmscut <= 0
        stdxyblim = max(stdxyb) + 1 ;
    else
        stdxyblim = spotrmscut ;
    end

    if spintb > 0 && max(stdxyb) < stdxyblim
        botspot = botcent4 ;
        bgood = 1 ;
    end

%    Find left spot

    [leftcent1,stdxy,spint] = ppspot(calim,leftcent,[segwidth segheight],cutlev) ;
    leftwidth = 2.*leftcent1(1) ;
    [leftcent2,stdxy,spint] = ppspot(calim,leftcent1,leftwidth,cutlev) ;
    if spotsize > leftwidth
        warning('Using smaller area than specified to find left spot: area falls outside image') ;
        leftspotsize = leftwidth ;
    else
        leftspotsize = spotsize ;
    end
    [leftcent3,stdxy,spint] = ppspot(calim,leftcent2,leftspotsize,cutlev) ;
    [leftcent4,stdxyl,spintl] = ppspot(calim,leftcent3,leftspotsize,cutlev) ;

    if spotrmscut <= 0
        stdxyllim = max(stdxyl) + 1 ;
    else
        stdxyllim = spotrmscut ;
    end

    if spintl > 0 && max(stdxyl) < stdxyllim
        leftspot = leftcent4 ;
        lgood = 1 ;
    end

%    Find right spot

    [rightcent1,stdxy,spint] = ppspot(calim,rightcent,[segwidth segheight],cutlev) ;
    rightwidth = 2.*( imwidth - rightcent1(1) ) ;
    [rightcent2,stdxy,spint] = ppspot(calim,rightcent1,rightwidth,cutlev) ;
    if spotsize > rightwidth
        warning('Using smaller area than specified to find right spot: area falls outside image') ;
        rightspotsize = rightwidth ;
    else
        rightspotsize = spotsize ;
    end
    [rightcent3,stdxy,spint] = ppspot(calim,rightcent2,rightspotsize,cutlev) ;
    [rightcent4,stdxyr,spintr] = ppspot(calim,rightcent3,rightspotsize,cutlev) ;

    if spotrmscut <= 0
        stdxyrlim = max(stdxyr) + 1 ;
    else
        stdxyrlim = spotrmscut ;
    end

    if spintr > 0 && max(stdxyr) < stdxyrlim
        rightspot = rightcent4 ;
        rgood = 1 ;
    end

%    Calculate spot spacing, image tilt and centre

    stspot = 'ttbbtl' ; enspot = 'lrlrbr' ;
    weighting = [1 1 1 1 2 2] ;
    lenf = sqrt(weighting) ;
    angmult = [-1 -1 -1 -1 1 -1] ;
    angoff = [pi./4 -pi./4 -pi./4 pi./4 0 0] ;
    alldists = 0 ; allangs = 0 ;
    weightdiv = 0 ;

    for i = 1:length(stspot)

        switch stspot(i)
            case 't'
                startspot = 'top' ;
            case 'b'
                startspot = 'bot' ;
            case 'l'
                startspot = 'left' ;
            case 'r'
                startspot = 'right' ;
        end

        switch enspot(i)
            case 't'
                finspot = 'top' ;
            case 'b'
                finspot = 'bot' ;
            case 'l'
                finspot = 'left' ;
            case 'r'
                finspot = 'right' ;
        end

        eval(['goodspots = ' stspot(i) 'good && ' enspot(i) 'good ;']) ;

        if goodspots

%    Calculate spot spacing for each combination

            eval([startspot finspot ' = ' startspot 'spot - ' finspot 'spot ;']) ;
            eval(['sfdist = ' startspot finspot ' ;']) ;
            eval([stspot(i) enspot(i) 'd = sqrt( sum( sfdist.^2 ) ) ;']) ;
            eval(['alldists = alldists + ' stspot(i) enspot(i) 'd.*lenf(i) ;']) ;
            weightdiv = weightdiv + weighting(i) ;

%    Calculate image tilt for each combination

            if strcmpi(stspot(i),'t') && strcmpi(enspot(i),'b')
                tba = ( atan(topbot(1)./topbot(2)) ) ;
            else
                eval([stspot(i) enspot(i) 'a = -( atan(' startspot finspot '(2)./' startspot finspot '(1)) + angoff(i) ) ;']) ;
            end
            eval(['allangs = allangs + (' stspot(i) enspot(i) 'a.*weighting(i) ) ;']) ;

%    Calculate centre prediction for each combination

            if strcmpi(stspot(i),'l') && strcmpi(enspot(i),'r')
                lrc = leftspot + (lrd./2).*[cos(lra) -sin(lra)] ;
            else
                if strcmpi(stspot(i),'b')
                    centmult = -1 ;
                else
                    centmult = 1 ;
                end
                eval([stspot(i) enspot(i) 'c = ' startspot 'spot + ' ... 
                    '(centmult.*(' stspot(i) enspot(i) 'd./( sqrt(2).*lenf(i) ) )).*[sin(' stspot(i) enspot(i) 'a) cos(' stspot(i) enspot(i) 'a)] ;']) ;
            end

            eval(['midpred(1) = midpred(1) + ' stspot(i) enspot(i) 'c(1).*weighting(i) ;']) ;
            eval(['midpred(2) = midpred(2) + ' stspot(i) enspot(i) 'c(2).*weighting(i) ;']) ;

        end

    end

%    Calculate spot spacing

    meandist = sqrt(2).*alldists./weightdiv ;
    onemil = meandist./spotdist ;
    sp2sp = onemil.*holepitch ;

%    Calculate image tilt

    imtilt = allangs./weightdiv ;

%    Find centre position

    if tgood && bgood && lgood && rgood

        m1 = topbot(2)./topbot(1) ;
        m2 = leftright(2)./leftright(1) ;
        A = [1 -m1; 1 -m2] ;
        c1 = topspot(2) - (topspot(1).*m1) ;
        c2 = leftspot(2) - (leftspot(1).*m2) ;
        c = [c1; c2] ;
        midpos = fliplr( (A\c)' ) ;

        if isinf(c1)
            midpos(1) = topspot(1) ;
        end

    else

        midpos = midpred./weightdiv ;

    end

%    Plot calibration image if so desired

    if doplot

        if makeplot
            figh = figure ;
        else
            figh = calfigh ;
        end

        figure(figh), imagesc(calim, [calmin calmax]), axis image ;
        colormap(gray) ;
        hold on ;
%        plot(topcent(1),topcent(2),'xy') ;
%        plot(botcent(1),botcent(2),'xy') ;
%        plot(leftcent(1),leftcent(2),'xy') ;
%        plot(rightcent(1),rightcent(2),'xy') ;
        if tgood
            topplot = topspot ;
            plot(topspot(1),topspot(2),'*r') ;
        else
            topplot = midpos - ( meandist./2 ).*[sin(imtilt) cos(imtilt)] ;
            plot(topplot(1),topplot(2),'*c') ;
        end
        if bgood
            botplot = botspot ;
            plot(botspot(1),botspot(2),'*r') ;
        else
            botplot = midpos + ( meandist./2 ).*[sin(imtilt) cos(imtilt)] ;
            plot(botplot(1),botplot(2),'*c') ;
        end
        if lgood
            leftplot = leftspot ;
            plot(leftspot(1),leftspot(2),'*r') ;
        else
            leftplot = midpos - ( meandist./2 ).*[cos(imtilt) -sin(imtilt)] ;
            plot(leftplot(1),leftplot(2),'*c') ;
        end
        if rgood
            rightplot = rightspot ;
            plot(rightspot(1),rightspot(2),'*r') ;
        else
            rightplot = midpos + ( meandist./2 ).*[cos(imtilt) -sin(imtilt)] ;
            plot(rightplot(1),rightplot(2),'*c') ;
        end
        plot([topplot(1),leftplot(1)], [topplot(2),leftplot(2)],'-r') ;
        plot([topplot(1),rightplot(1)], [topplot(2),rightplot(2)],'-r') ;
        plot([botplot(1),leftplot(1)], [botplot(2),leftplot(2)],'-r') ;
        plot([botplot(1),rightplot(1)], [botplot(2),rightplot(2)],'-r') ;
        plot([topplot(1),botplot(1)], [topplot(2),botplot(2)],'-r') ;
        plot([leftplot(1),rightplot(1)], [leftplot(2),rightplot(2)],'-r') ;
        plot(midpos(1),midpos(2),'*c') ;
        hold off ;

    end

%    Change sign of tilt angle to match pepperpot analysis code

    imtilt = -imtilt ;

%    Display calibration data

    if nargout < 1
        disp(' ') ;
        disp(['Centre position of image = [' num2str(midpos(1)) ' ' num2str(midpos(2)) '] pixels.']) ;
        disp(' ') ;
        disp(['Spot-to-spot spacing = ' num2str(sp2sp) ' pixels.']) ;
        disp(' ') ;
        disp(['Angle of image tilt = ' num2str(imtilt) ' radians.']) ;
        disp(' ') ;
    end

    tout = topspot ; bout = botspot ;
    lout = leftspot ; rout = rightspot ;

    return

%    Original calculations done by hand

    if 0

        topleft = topspot - leftspot ;
        topright = topspot - rightspot ;
        botleft = botspot - leftspot ;
        botright = botspot - rightspot ;
        topbot = topspot - botspot ;
        leftright = leftspot - rightspot ;

        tld = sqrt( sum( topleft.^2 ) ) ;
        trd = sqrt( sum( topright.^2 ) ) ;
        bld = sqrt( sum( botleft.^2 ) ) ;
        brd = sqrt( sum( botright.^2 ) ) ;
        tbd = sqrt( sum( topbot.^2 ) ) ;
        lrd = sqrt( sum( leftright.^2 ) ) ;

        meandist = sqrt(2).*( tld + trd + bld + brd + ( sqrt(2).*(tbd + lrd) ) )./8 ;

%    Calculate image tilt

        tla = -( atan(topleft(2)./topleft(1)) + pi./4 ) ;
        tra = -( atan(topright(2)./topright(1)) - pi./4 ) ;
        bla = -( atan(botleft(2)./botleft(1)) - pi./4 ) ;
        bra = -( atan(botright(2)./botright(1)) + pi./4 ) ;
        tba = ( atan(topbot(1)./topbot(2)) ) ;
        lra = -( atan(leftright(2)./leftright(1)) ) ;
        imtilt = ( tla + tra + bla + bra + 2.*( tba + lra ) )./8 ;

%    Find centre position

        tlc = topspot + (tld./sqrt(2)).*[sin(tla) cos(tla)] ;
        trc = topspot + (trd./sqrt(2)).*[sin(tra) cos(tra)] ;
        blc = botspot - (tld./sqrt(2)).*[sin(bla) cos(bla)] ;
        brc = botspot - (trd./sqrt(2)).*[sin(bra) cos(bra)] ;
        tbc = topspot + (tbd./2).*[sin(tba) cos(tba)] ;
        lrc = leftspot + (lrd./2).*[cos(lra) -sin(lra)] ;

    end




