% PPCALIB
%
%    PPCALIB - Calibrate pepperpot images

    clear all ; close all ;

    [ppim,stdm,pphist] = pepperplotfun('E:\pepperpot\04_04_06\cam_back_screw-shims.tif',0,[2 0]) ;

%    figure, imagesc(ppim(:,:),[1 16384]), axis square, hold on ;
%    plot(999.5,1007.5,'xr') ;   % Mark central point
%    axis tight ;
%    axis([970 1026 980 1034]) ;
%    hold off ;

    ppbw = edge(ppim,'canny') ;
    figure, imshow(ppbw) ;

%    cent = [999.5 1007.5] ;
    cent = [1000 1008] ;
    rull = [113 1884] ;
    rulr = [1862 1895] ;

%    figure, imagesc(ppim(:,:),[1 16384]), axis square, hold on ;
%    colormap(gray) ;
%    plot([rull(1),rulr(1)], [rull(2),rulr(2)],'-xr') ;
%    plot([rull(1),cent(1)], [rull(2),cent(2)],'-xr') ;
%    plot([cent(1),rulr(1)], [cent(2),rulr(2)],'-xr') ;
%    hold off ;

    cldiff = abs( cent - rull ) ;   % Difference between left ruler end and central point
    crdiff = abs( rulr - cent ) ;   % Difference between right ruler end and central point
    ruldiff = abs( rulr - rull ) ;  % Difference between end points of ruler
    rulgrad = abs( ruldiff(2) ./ ruldiff(1) ) ;     % Gradient or ruler
    ruloff = rull(2) - ( rull(1).*ruldiff(2) ) ;    % Offset of ruler
    rulmid = rull + ruldiff./2 ;                    % Midpoint of ruler
    mcdiff = abs( rulmid - cent ) ;     % Difference between ruler midpoint and centre

    rulen = sqrt( sum(ruldiff.^2) ) ;   % Length of ruler
    cllen = sqrt( sum(cldiff.^2) ) ;    % Distance between central point and left ruler edge
    crlen = sqrt( sum(crdiff.^2) ) ;    % Distance between central point and right ruler edge
    mclen = sqrt( sum(mcdiff.^2) ) ;    % Distance between central point and ruler midpoint
    ruang = atan( ruldiff(2)./ruldiff(1) ) ;    % Angle of ruler (and therefore image)
    clang = atan( cldiff(2)./cldiff(1) ) ;      % Angle of Centre-to-ruler-left line
    crang = atan( crdiff(2)./crdiff(1) ) ;      % Angle of Centre-to-ruler-right line
    mcang = atan( mcdiff(2)./mcdiff(1) ) ;      % Angle of midpoint-to-centre line

    lang = clang + ruang ;
    rang = crang - ruang ;
    cang = (pi - lang - rang ) ;
    mang = mcang + ruang ;

%    Refind centre position
%    Angle rule - right corner point

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

    coffr = [ ( rlen .* cos(brang) ), ( rlen .* sin(brang) ) ] ;
    cposr = rulr - coffr ;

%    Angle rule - left corner point

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

    coffl = [ -( llen .* cos(blang) ), ( llen .* sin(blang) ) ] ;
    cposl = rull - coffl ;

%    Replot mid position from both measurements

    figure, imagesc(ppim, [1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot([rull(1),rulr(1)], [rull(2),rulr(2)],'-xr') ;
    plot([rull(1),cent(1)], [rull(2),cent(2)],'-xr') ;
    plot([cent(1),rulr(1)], [cent(2),rulr(2)],'-xr') ;
    plot(cposr(1),cposr(2),'*c') ;
    plot(cposl(1),cposl(2),'+y') ;
    plot(rulmid(1),rulmid(2),'xc') ;

%    Midline rule

    midp = rulmid ;  mpang = mang - ruang ;
    mcoff = [ -( mclen .* cos(mpang) ), ( mclen .* sin(mpang) ) ] ;
    cposm = midp - mcoff ;
    plot(midp(1),midp(2),'xc') ;
    plot([midp(1),cposm(1)], [midp(2),cposm(2)],'-xc') ;
    hold off ;

    rhspot = [1764 1011] ;
    lhspot = [238 1005] ;
    mltilt = atan2( (rhspot(2) - lhspot(2)),(rhspot(1) - lhspot(1)) ) ;
    rumlang = mltilt - ruang ;

    figure, imagesc(ppim, [1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot([rull(1),rulr(1)], [rull(2),rulr(2)],'-xr') ;
    plot([rull(1),cent(1)], [rull(2),cent(2)],'-xr') ;
    plot([cent(1),rulr(1)], [cent(2),rulr(2)],'-xr') ;
    plot([lhspot(1),rhspot(1)], [lhspot(2),rhspot(2)],'-xr') ;
    plot(cposr(1),cposr(2),'*c') ;
    plot([rulmid(1),cposm(1)], [rulmid(2),cposm(2)],'-xc') ;
    hold off ;

%    Use ppcalimdat function to measure data again

    [outdat,outnames] = ppcalimdat(ppim,cent,rull,rulr,lhspot,rhspot) ;

%    Test position finder with new image

    [ppim2,stdm2,pphist2] = pepperplotfun('E:\pepperpot\04_04_06\cam_in_screw-shims.tif',0,[2 0]) ;

    figure, imagesc(ppim2, [1 16384]), axis square ;
    colormap(gray) ;

    ppbw2 = edge(ppim2,'canny') ;
    figure, imagesc(ppbw2), axis square ;
    colormap(gray) ;

    cent2 = [1134 1067.5] ;
    rull2 = [251 1943] ;
    rulr2 = [1998 1953] ;

    figure, imshow(ppim2, [1 16384]), hold on ;
    plot([rull2(1),rulr2(1)], [rull2(2),rulr2(2)],'-xr') ;
    plot([rull2(1),cent2(1)], [rull2(2),cent2(2)],'-xr') ;
    plot([cent2(1),rulr2(1)], [cent2(2),rulr2(2)],'-xr') ;
    hold off ;

    ruldiff2 = abs( rulr2 - rull2 ) ;  % Difference between end points of ruler
    rulen2 = sqrt( sum(ruldiff2.^2) ) ;   % Length of ruler
    ruang2 = atan( ruldiff2(2)./ruldiff2(1) ) ;    % Angle of ruler (and therefore image)
    rulmid2 = rull2 + (ruldiff2./2) ;   % Centre point of ruler

%    Refind centre position
%    Angle rule - right corner point

    rlen2 = rulen2 .* ( sin(lang) ./ sin(cang) ) ;
    brang2 = rang + ruang2 ;

    coffr2 = [ ( rlen2 .* cos(brang2) ), ( rlen2 .* sin(brang2) ) ] ;
    cposr2 = round(rulr2 - coffr2) ;

%    Angle rule - left corner point

    llen2 = rulen2 .* ( sin(rang) ./ sin(cang) ) ;
    blang2 = lang - ruang2 ;

    coffl2 = [ -( llen2 .* cos(blang2) ), ( llen2 .* sin(blang2) ) ] ;
    cposl2 = round(rull2 - coffl2) ;

%    Midline rule

    mpang2 = mang - ruang2 ;
    mcoff2 = [ -( mclen .* cos(mpang2) ), ( mclen .* sin(mpang2) ) ] ;
    cposm2 = rulmid2 - mcoff2 ;

%    Replot mid position from both measurements

    figure, imshow(ppim2, [1 16384]), hold on ;
    plot([rull2(1),rulr2(1)], [rull2(2),rulr2(2)],'-xr') ;
    plot([rull2(1),cent2(1)], [rull2(2),cent2(2)],'-xr') ;
    plot([cent2(1),rulr2(1)], [cent2(2),rulr2(2)],'-xr') ;
    plot(cposr2(1),cposr2(2),'*c') ;
    plot(cposl2(1),cposl2(2),'+y') ;
    plot(rulmid2(1),rulmid2(2),'xc') ;
    plot([rulmid2(1),cposm2(1)], [rulmid2(2),cposm2(2)],'-xc') ;
    hold off ;

    rhspot2 = [1897 1070] ;
    lhspot2 = [374 1066] ;
    mltilt2 = atan2( (rhspot2(2) - lhspot2(2)),(rhspot2(1) - lhspot2(1)) ) ;

    figure, imagesc(ppim2, [1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot([rull2(1),rulr2(1)], [rull2(2),rulr2(2)],'-xr') ;
    plot([rull2(1),cent2(1)], [rull2(2),cent2(2)],'-xr') ;
    plot([cent2(1),rulr2(1)], [cent2(2),rulr2(2)],'-xr') ;
    plot([lhspot2(1),rhspot2(1)], [lhspot2(2),rhspot2(2)],'-xr') ;
    plot(cposr2(1),cposr2(2),'*c') ;
    plot([rulmid2(1),cposm2(1)], [rulmid2(2),cposm2(2)],'-xc') ;
    hold off ;

    rumlang2 = mltilt2 - ruang2 ;

%    Use ppcalimdat function to measure data again

    [outdat2] = ppcalimdat(ppim2,cent2,rull2,rulr2,lhspot2,rhspot2) ;

%    Test position finder with new image

    [ppim3,stdm3,pphist3] = pepperplotfun('E:\pepperpot\03_04_06\cam_back_shimmed.tif',0,[2 0]) ;

    cent3 = [994 1035] ;
    rull3 = [99 1912] ;
    rulr3 = [1856 1931] ;
    lhspot3 = [229 1029] ;
    rhspot3 = [1762 1043] ;

%    Use ppcalimdat function to measure data

    [outdat3] = ppcalimdat(ppim3,cent3,rull3,rulr3,lhspot3,rhspot3) ;

%    Test position finder with new image

    [ppim4,stdm4,pphist4] = pepperplotfun('E:\pepperpot\03_04_06\cam_in_shimmed.tif',0,[2 0]) ;

    cent4 = [1021 1140] ;
    rull4 = [128 2017] ;
    rulr4 = [1884 2034] ;
    lhspot4 = [257 1135] ;
    rhspot4 = [1788 1146] ;

%    Use ppcalimdat function to measure data

    [outdat4] = ppcalimdat(ppim4,cent4,rull4,rulr4,lhspot4,rhspot4) ;

    rulenm = mean([outdat{6} outdat2{6} outdat3{6} outdat4{6}]) ;
    langm = mean([outdat{8} outdat2{8} outdat3{8} outdat4{8}]) ;
    rangm = mean([outdat{9} outdat2{9} outdat3{9} outdat4{9}]) ;
    cangm = mean([outdat{10} outdat2{10} outdat3{10} outdat4{10}]) ;
    rumlangm = mean([outdat{12} outdat2{12} outdat3{12} outdat4{12}]) ;
    mllenm = mean([outdat{13} outdat2{13} outdat3{13} outdat4{13}]) ;
    mlrulratio = mllenm./rulenm ;

    rupix = rulenm./80 ;    % Number of pixels per mm for ruler image
    sppix = mllenm./72 ;   % Number of pixels per mm for spot image
    pixscale = rupix./sppix ;
    phosd = 13 ;    % Distance to phosphor from tungsten screen
    ruld = 19 ;     % Distance to rulers from tungsten screen
    fracdist = phosd./ruld ;    % Fractional distance to phosphor screen
    gradscale = ( pixscale - 1 ) ;
    perscale = gradscale.*fracdist + 1 ;
    invperscale = 1./perscale ;

%    save('E:\pepperpot\cal_data\pixcal.mat','rupix','perscale','invperscale','mlrulratio') ;

    [midpos1,imtilt1,rulen1a] = pptiltfind(ppim,rull,rulr) ;
    [midpos2,imtilt2,rulen2a] = pptiltfind(ppim2,rull2,rulr2) ;
    [midpos3,imtilt3,rulen3a] = pptiltfind(ppim3,rull3,rulr3) ;
    [midpos4,imtilt4,rulen4a] = pptiltfind(ppim4,rull4,rulr4) ;

%    load E:\pepperpot\cal_data\ppcaldat2.mat ;
%    ppcalmat(:,:,1) = ppcalmat(:,:,1).*perscale.*mlrulratio ;
%    save('E:\pepperpot\cal_data\ppcaldat_good.mat','ppcalmat') ;
    load E:\pepperpot\cal_data\ppcaldat_good.mat ;

    [ppimreco] = ppspotreco(ppcalmat,midpos1,rulen1a,imtilt1,invperscale) ;

    figure, imagesc(ppim,[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(midpos1(1),midpos1(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppimreco(i,j,1),ppimreco(i,j,2),'xc') ;
        end
    end
    hold off

    [ppim2reco] = ppspotreco(ppcalmat,midpos2,rulen2a,imtilt2,invperscale) ;

    figure, imagesc(ppim2,[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(midpos2(1),midpos2(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppim2reco(i,j,1),ppim2reco(i,j,2),'xc') ;
        end
    end
    hold off

    [ppim3reco] = ppspotreco(ppcalmat,midpos3,rulen3a,imtilt3,invperscale) ;

    figure, imagesc(ppim3,[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(midpos3(1),midpos3(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppim3reco(i,j,1),ppim3reco(i,j,2),'xc') ;
        end
    end
    hold off

    [ppim4reco] = ppspotreco(ppcalmat,midpos4,rulen4a,imtilt4,invperscale) ;

    figure, imagesc(ppim4,[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(midpos4(1),midpos4(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppim4reco(i,j,1),ppim4reco(i,j,2),'xc') ;
        end
    end
    hold off

    return


