function [midpos,sp2sp,imtilt] = ppcalfindmm(immat,tleft,tright,bleft,bright,mmsize,plt,calfigh)
% Function PPCALFINDMM
%
%    [midpos,sp2sp,imtilt] = ppcalfindmm(immat,tleft,tright,bleft,bright,mmsize,plt,calfigh)
%
%    PPCALFINDMM is used to find the rotation, centre and scaling of
%    Pepperpot data images for a particular data set, using the position of
%    4 calibration markings on the copper mounting block, each spaced 9mm,
%    63mm or 126mm apart, on a calibration image.
%
%    PPCALFINDMM(IMMAT) loads the image matrix or the image file IMMAT and
%    displays it on screen: the user can then find the pixel positions of
%    the top-left, top-right, bottom-left and bottom-right calibration
%    markings and the spacing of those markings (9, 63 or 126mm) and input
%    this information at the prompt in the form of a matrix:
%
%       > Please enter top-left marker coordinates:   [950 950]
%
%      The centre of the image, the rotation of the image and the
%      spot-to-spot spacing are then calculated and displayed on screen.
%
%    If IMMAT is a string, PPCALFINDMM calls PEPPERPLOTFUN to load the
%    image data; otherwise it is assumed to be a matlab variable.
%
%    [MIDPOS] = PPCALFINDMM(IMMAT) outputs the centre position of the
%    calibration image to the variable MIDPOS.
%
%    [MIDPOS,SP2SP] = PPCALFINDMM(IMMAT) also outputs the spot-to-spot
%    spacing of the mesh, in pixels, to the variable SP2SP.
%
%    [MIDPOS,SP2SP,IMTILT] = PPCALFINDMM(IMMAT) also outputs the tilt of the
%    image in radians to the variable IMTILT.
%
%    PPCALFINDMM(IMMAT,[TLEFT],[TRIGHT],[BLEFT],[BRIGHT],MMSIZE) makes it
%    possible to specify the four pixel locations of the top-left,
%    top-right, bottom-left and bottom-right calibration markings, where: 
%
%        TLEFT = [TLX TLY] ;    TRIGHT = [TRX TRY] ;
%        BLEFT = [BLX BLY] ;    BRIGHT = [TRX TRY] ;
%
%      ... and the marker spacing, MMSIZE, in millimetres: this must either
%    be 9, 63 or 126.  The coordinate positions of ALL FOUR MARKINGS must
%    be specified, as well as the marker spacing, otherwise the function
%    will not run.
%
%    PPCALFINDMM(IMMAT,[TLEFT],[TRIGHT],[BLEFT],[BRIGHT],MMSIZE,PLT) controls
%    whether the background image is plotted: it PLT = 1, the data and fit
%    are plotted.  The default value is 1 (data is plotted).
%
%    PPCALFINDMM(IMMAT,[TLEFT],[TRIGHT],[BLEFT],[BRIGHT],MMSIZE,PLT,CALFIGH) 
%    plots the data in figure(CALFIGH) to allow GUI plotting.
%
%    The function requires the following input variables:
%
%    IMMAT - The data/image file with the image information, including
%              the file extension and the path (absolute or relative) -
%              must be a string to load data or a variable in the Matlab
%              workspace that already contains the image data.
%
%    It is possible to specify IMMAT as an empty matrix to carry out the
%    calibration calculations without actually loading the image ie.
%
%    [MIDPOS,SP2SP,IMTILT] = PPCALFINDMM([],[TLEFT],[TRIGHT],[BLEFT],[BRIGHT],MMSIZE,PLT)
%
%    In this case, the PLT variable is always set to zero.
%

    if nargin < 1
        error('Need at least 1 input variable: [midpos,sp2sp,imtilt] = ppcalfindmm(immat,tleft,tright,bleft,bright,mmsize,plt)') ;
    elseif nargin > 8
        error('Can only specify 8 input variables: [midpos,sp2sp,imtilt] = ppcalfindmm(immat,tleft,tright,bleft,bright,mmsize,plt,calfigh)') ;
    elseif ( nargin == 2 ) || ( nargin == 3 ) || ( nargin == 4 )
        error('Must specify coordinates for all 4 markings: ppcalfindmm(immat,tleft,tright,bleft,bright,mmsize,plt)') ;
    elseif nargout > 3
        error('Can only specify 3 output variables: [midpos,sp2sp,imtilt] = ppcalfindmm(immat,tleft,tright,bleft,bright,mmsize,plt)') ;
    elseif nargin > 4
        if ( length(tleft) ~= 2 ) || ( length(tright) ~= 2 ) || ( length(bleft) ~= 2 ) || ( length(bright) ~= 2 )
            error('Ruler coordinates must contain X and Y pixels: tleft = [tlx tly]') ;
        end
    elseif ( isempty(immat) ) && ( nargin < 5 )
        error('Must specify all calibration locations if using empty image matrix: ppcalfindmm([],tleft,tright,bleft,bright)') ;
    end

    if nargin < 7
        plt = 1 ;
    end

    if nargin < 6
        mmsize = 9 ;
    end

    if ischar(immat)

%        if plt == 1
%            funplt = 2 ;
%        else
            funplt = 0 ;
%        end

        [ppim] = pepperplotfun(immat,0,[funplt 0]) ;
        immat = ppim ;

    end

    if ~isempty(immat)
        immin = min(min(immat)) ; 
        immax = max(max(immat)) ;
    end

%    Ask for positions of ruler ends if not provided and plot data

    if nargin < 2

        gsim = figure ; 

        if ~isempty(immat)
            imagesc(immat, [immin immax]), axis image ;
            colormap(gray) ;
        end

%        ppbw = edge(immat,'canny') ;
%        bwim = figure ; imagesc(ppbw), axis square ;
%        colormap(gray) ;

        disp(' ') ;
        tleft = input('Please enter top-left marking coordinates:  ') ;
        disp(' ') ;
        tright = input('Please enter top-right marking coordinates:  ') ;
        disp(' ') ;
        bleft = input('Please enter bottom-left marking coordinates:  ') ;
        disp(' ') ;
        bright = input('Please enter bottom-right marking coordinates:  ') ;
        disp(' ') ;
        mmsize = input('Please enter spacing of marking in mm:  ') ;
        disp(' ') ;

        figure(gsim), hold on ;
        plot([tleft(1),tright(1)], [tleft(2),tright(2)],'-xr') ;
        plot([tleft(1),bleft(1)], [tleft(2),bleft(2)],'-xr') ;
        plot([tright(1),bright(1)], [tright(2),bright(2)],'-xr') ;
        plot([bleft(1),bright(1)], [bleft(2),bright(2)],'-xr') ;
        plot([tleft(1),bright(1)], [tleft(2),bright(2)],'-xr') ;
        plot([bleft(1),tright(1)], [bleft(2),tright(2)],'-xr') ;
        hold off ;

    end

    goodsizes = [9 63 126 132] ;

    if ~sum( mmsize == goodsizes )
        error('Marker spacing MMSIZE must be 9, 63 or 126 mm') ;
    end

    markspac = mmsize./3 ;

    tltr = tleft - tright ;
    tlbl = tleft - bleft ;
    trbr = tright - bright ;
    blbr = bleft - bright ;
    tlbr = tleft - bright ;
    trbl = tright - bleft ;

    tltrd = sqrt( sum( tltr.^2 ) ) ;
    tlbld = sqrt( sum( tlbl.^2 ) ) ;
    trbrd = sqrt( sum( trbr.^2 ) ) ;
    blbrd = sqrt( sum( blbr.^2 ) ) ;
    tlbrd = sqrt( sum( tlbr.^2 ) ) ;
    trbld = sqrt( sum( trbl.^2 ) ) ;
    rt2 = sqrt(2) ;

    meandist = ( tltrd + tlbld + trbrd + blbrd )./4 ;
    sp2sp = meandist./markspac ;

%    Find centre position

    m1 = trbl(2)./trbl(1) ;
    m2 = tlbr(2)./tlbr(1) ;
    A = [1 -m1; 1 -m2] ;
    c1 = tright(2) - (tright(1).*m1) ;
    c2 = tleft(2) - (tleft(1).*m2) ;
    c = [c1; c2] ;
    midpos = fliplr( (A\c)' ) ;

%    Calculate image tilt

    tltra = atan(tltr(2)./tltr(1)) ;
    blbra = atan(blbr(2)./blbr(1)) ;
    tlbla = -( atan(tlbl(1)./tlbl(2)) ) ;
    trbra = -( atan(trbr(1)./trbr(2)) ) ;
    imtilt = ( tltra + blbra + tlbla + trbra )./4 ;

%    Plot mid position

    if ( nargin > 7 ) || ( ( plt == 1 ) && ( ~isempty(immat) ) )

        if nargin > 7
            figure(calfigh) ; hold on ;
        else
            figure, imagesc(immat, [immin immax]), axis image, hold on ;
            colormap(gray) ;
        end
        plot([tleft(1),tright(1)], [tleft(2),tright(2)],'-xr') ;
        plot([tleft(1),bleft(1)], [tleft(2),bleft(2)],'-xr') ;
        plot([tright(1),bright(1)], [tright(2),bright(2)],'-xr') ;
        plot([bleft(1),bright(1)], [bleft(2),bright(2)],'-xr') ;
        plot([tleft(1),bright(1)], [tleft(2),bright(2)],'-xr') ;
        plot([bleft(1),tright(1)], [bleft(2),tright(2)],'-xr') ;
        plot(midpos(1),midpos(2),'*c') ;
        hold off ;

    end

%    Display calibration data

%    midpos = round(cposr) ;
%    imtilt = ( ruang + rumlang ) ;

    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(' ') ;

    return








