function [midpos,imtilt,rulen] = pptiltfind(immat,lpos,rpos,plt)
% Function PPTILTFIND
%
%    [midpos,imtilt,rulen] = pptiltfind(immat,lpos,rpos,plt)
%
%    PPTILTFIND is used to find the rotation and centre of Pepperpot data
%    images for a particular data set, using the position of the horizontal
%    ruler on a calibration image.
%
%    PPTILTFIND(IMMAT) loads the image matrix IMMAT and displays it on
%    screen: the user can then find the left and right hand edges of the
%    horizontal ruler by hand and input this information at the prompt in
%    the form of a matrix:
%
%       > Please enter left-hand ruler coordinates:   [113 1884]
%
%      The centre of the image and the rotation of the image are then
%    calculated and displayed on screen.
%
%    The image data should be preloaded using the PEPPERPLOTFUN function.
%
%    [MIDPOS] = PPTILTFIND(IMMAT) outputs the centre position of the
%    calibration image to the variable MIDPOS.
%
%    [MIDPOS,IMTILT] = PPTILTFIND(IMMAT) also outputs the tilt of the
%    image in radians to the variable IMTILT.
%
%    PPTILTFIND(IMMAT,[LPOS],[RPOS]) makes it possible to specify the
%    pixel locations of the left- and right-hand ends of the ruler, where:
%
%        LPOS = [LX LY] ;    RPOS = [RX RY] ;
%
%      The coordinate positions of BOTH ends of the ruler must be specified
%    otherwise the function will not run.
%
%    PPTILTFIND(IMMAT,[LPOS],[RPOS],PLT) controls whether the background
%    image are plotted: it PLT = 1, the data and fit are plotted.  The
%    default value is 1 (data is plotted).
%
%    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.
%

    if nargin < 1
        error('Need at least 1 input variable: [midpos,imtilt,rulen] = pptiltfind(immat,lpos,rpos,plt)') ;
    elseif nargin > 4
        error('Can only specify 4 input variables: [midpos,imtilt,rulen] = pptiltfind(immat,lpos,rpos,plt)') ;
    elseif nargin == 2
        error('Must specify coordinates for both ends of ruler: [midpos,imtilt,rulen] = pptiltfind(immat,lpos,rpos,plt)') ;
    elseif nargout > 3
        error('Can only specify 3 output variables: [midpos,imtilt,rulen] = pptiltfind(immat,lpos,rpos,plt)') ;
    elseif nargin == 3
        if ( length(lpos) ~= 2 ) | ( length(rpos) ~= 2 )
            error('Ruler coordinates must contain X and Y pixels: lpos = [lx ly]') ;
        end
    end

    if nargin < 4
        plt = 1 ;
    end

%    Data from original calibration images

    lang = -0.786099737270937 ;         % Angle between left-hand triangle edge and ruler
    rang = 0.793161570054275 ;          % Angle between right-hand triangle edge and ruler
    cang = 1.56233134626458 ;           % Angle between left- and right-hand triangle edges
    rumlang = -0.00240810298182418 ;    % Angle between ruler and image midline

    immin = min(min(immat)) ; immax = max(max(immat)) ;

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

    if nargin < 2

        gsim = figure ; imagesc(immat, [immin immax]), axis image ;
        colormap(gray) ;

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

        disp(' ') ;
        lpos = input('Please enter left-hand ruler coordinates:  ') ;
        disp(' ') ;
        rpos = input('Please enter right-hand ruler coordinates:  ') ;
        disp(' ') ;

        figure(gsim), hold on ;
        plot([lpos(1),rpos(1)], [lpos(2),rpos(2)],'-xr') ;
        hold off ;

    end

    ruldiff = abs( rpos - lpos ) ;  % Difference between end points of ruler
    rulen = sqrt( sum(ruldiff.^2) ) ;   % Length of ruler
    ruang = atan( ruldiff(2)./ruldiff(1) ) ;    % Angle of ruler (and therefore image)
    rulmid = lpos + (ruldiff./2) ;   % Centre point of ruler

%    Find centre position
%    Angle rule - right corner point

    rlen = rulen .* ( sin(lang) ./ sin(cang) ) ;
    brang = rang + ruang ;

    coffr = [ ( rlen .* cos(brang) ), ( rlen .* sin(brang) ) ] ;
    cposr = round(rpos + coffr) ;

%    Angle rule - left corner point

    llen = rulen .* ( sin(rang) ./ sin(cang) ) ;
    blang = lang + ruang ;

    coffl = [ ( llen .* cos(blang) ), ( llen .* sin(blang) ) ] ;
    cposl = round(lpos + coffl) ;

%    Plot mid position from both measurements

    if plt == 1

        figure, imagesc(immat, [immin immax]), axis image, hold on ;
        colormap(gray) ;
        plot([lpos(1),rpos(1)], [lpos(2),rpos(2)],'-xr') ;
        plot([lpos(1),cposl(1)], [lpos(2),cposl(2)],'-xr') ;
        plot([cposr(1),rpos(1)], [cposr(2),rpos(2)],'-xr') ;
        plot(cposr(1),cposr(2),'*c') ;
        plot(cposl(1),cposl(2),'+y') ;
        plot(rulmid(1),rulmid(2),'xc') ;
        hold off ;

    end

%    Display calibration data

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

    disp(' ') ;
    disp(['Centre position of image = ' num2str(midpos) ' pixels.']) ;
    disp(' ') ;
    disp(['Angle of image tilt = ' num2str(imtilt) ' radians.']) ;
    disp(' ') ;

    return








