% PPFINDHOLES_MK2MESH
%
%    PPFINDHOLES_MK2MESH - Create calibration matrix using MkII pepperpot
%    mesh image

    clear all ; close all ;

%    [ppim,stdm,pphist,ppminmax] = pepperplotfun('E:\pepperpot\41x41meshscan\ppmesh_scan1.tif',0,[2 0]) ;
%    ppmin = ppminmax(1) ; ppmax = ppminmax(2) ;
%    ppim = imread('E:\pepperpot\41x41meshscan\ppmesh_scan1.tif','tif') ;
%    ppim = ppim(251:14050, 351:14150) ;

%    ppimi = imread('E:\pepperpot\41x41meshscan\ppmesh_scan1_cut.tif','tif') ;
%    ppim = single(ppimi) ;

    load E:\pepperpot\41x41meshscan\ppmesh_scan1_dat.mat ppim ;
    ppmin = min(min(ppim)) ; ppmax = max(max(ppim)) ;

    pplim = 5000 ;
%    pplim2 = 500 ;
%    pplimlow = 130 ;
%    ppcut = ( ppim > pplim ) ;
%    ppcut2 = ( ppim > pplim2 ) ;
%    ppimcut = ( ppim - pplim ).*ppcut ;

%    figure, imagesc(ppcut), axis image, hold on ;
%    colormap(gray) ;
%    hold off ;

    tlb = [256 266] ;
    trb = [13504 295] ;
    blb = [253 13481] ;
    brb = [13480 13503] ;
    rb = [13505 6892] ;
    lb = [245 6866] ;
    tb = [6910 280] ;
    bb = [6902 13492] ;

    cent = [6906 6884] ;

    tl = [1205 1214] ;
    tr = [12563 1238] ;
    bl = [1203 12539] ;
    br = [12561 12560] ;
    r = [12570 6895] ;
    l = [1198 6873] ;
    t = [6908 1224] ;
    b = [6903 12549] ;

    if 0
        figure, imagesc(ppimi(:,:),[ppmin ppmax]), axis image, hold on ;
        colormap(gray) ;
        plot(cent(1),cent(2),'xr') ;    % Mark central point
        plot(tlb(1),tlb(2),'xy') ;        % Mark upper left support hole
        plot(trb(1),trb(2),'xy') ;        % Mark upper right support hole
        plot(blb(1),blb(2),'xy') ;        % Mark lower left support hole
        plot(brb(1),brb(2),'xy') ;        % Mark lower right support hole
        plot(lb(1),lb(2),'xy') ;          % Mark left support hole
        plot(tb(1),tb(2),'xy') ;          % Mark top support hole
        plot(bb(1),bb(2),'xy') ;          % Mark bottom support hole
        plot(rb(1),rb(2),'xy') ;          % Mark right support hole
        plot([tl(1),tr(1)], [tl(2),tr(2)],'-xr') ;
        plot([tl(1),bl(1)], [tl(2),bl(2)],'-xr') ;
        plot([tr(1),br(1)], [tr(2),br(2)],'-xr') ;
        plot([bl(1),br(1)], [bl(2),br(2)],'-xr') ;
        plot([tl(1),br(1)], [tl(2),br(2)],'-xr') ;
        plot([tr(1),bl(1)], [tr(2),bl(2)],'-xr') ;
        plot([t(1),b(1)], [t(2),b(2)],'-xr') ;
        plot([l(1),r(1)], [l(2),r(2)],'-xr') ;
        hold off ;
    end

%    axis([1190 1210 12530 12550]) ;
%    axis([12550 12570 12550 12570]) ;
%    axis([1195 1215 1203 1223]) ;
%    axis([12550 12570 1230 1250]) ;
%    axis([6850 6950 6840 6940]) ;

    spotwidth = 41 ;
    spotr = (spotwidth - 1)./2 ;
    intspotw = 21 ;

    [cent12,cent12std] = ppspot(ppim,cent,spotwidth,pplim) ;
    [br12,br12std] = ppspot(ppim,br,spotwidth,pplim) ;
    [bl12,bl12std] = ppspot(ppim,bl,spotwidth,pplim) ;
    [tr12,tr12std] = ppspot(ppim,tr,spotwidth,pplim) ;
    [tl12,tl12std] = ppspot(ppim,tl,spotwidth,pplim) ;
    [b12,b12std] = ppspot(ppim,b,spotwidth,pplim) ;
    [l12,l12std] = ppspot(ppim,l,spotwidth,pplim) ;
    [r12,r12std] = ppspot(ppim,r,spotwidth,pplim) ;
    [t12,t12std] = ppspot(ppim,t,spotwidth,pplim) ;

    if 0
        figure, imagesc(ppimi(:,:),[ppmin ppmax]), axis image, hold on ;
        colormap(gray) ;
        plot(cent(1),cent(2),'xr') ;    % Mark central point
        plot(tl(1),tl(2),'xr') ;        % Mark upper left point
        plot(tr(1),tr(2),'xr') ;        % Mark upper right point
        plot(bl(1),bl(2),'xr') ;        % Mark lower left point
        plot(br(1),br(2),'xr') ;        % Mark lower right point
        plot(l(1),l(2),'xr') ;          % Mark left point
        plot(t(1),t(2),'xr') ;          % Mark top point
        plot(b(1),b(2),'xr') ;          % Mark bottom point
        plot(r(1),r(2),'xr') ;          % Mark right point
        plot(cent12(1),cent12(2),'xc') ;    % Mark new central point
        plot(tl12(1),tl12(2),'xc') ;      % Mark new upper left point
        plot(tr12(1),tr12(2),'xc') ;      % Mark new upper right point
        plot(bl12(1),bl12(2),'xc') ;      % Mark new lower left point
        plot(br12(1),br12(2),'xc') ;      % Mark new lower right point
        plot(l12(1),l12(2),'xc') ;          % Mark new left point
        plot(t12(1),t12(2),'xc') ;          % Mark new top point
        plot(b12(1),b12(2),'xc') ;          % Mark new bottom point
        plot(r12(1),r12(2),'xc') ;          % Mark new right point
        hold off ;
    end

    if 0
        figure, imagesc(ppimi(:,:),[ppmin ppmax]), axis image, hold on ;
        colormap(gray) ;
        plot(cent(1),cent(2),'xr') ;    % Mark central point
        plot(cent12(1),cent12(2),'xc') ;    % Mark new central point
        plot([tl(1),tr(1)], [tl(2),tr(2)],'-xr') ;
        plot([tl(1),bl(1)], [tl(2),bl(2)],'-xr') ;
        plot([tr(1),br(1)], [tr(2),br(2)],'-xr') ;
        plot([bl(1),br(1)], [bl(2),br(2)],'-xr') ;
        plot([tl(1),br(1)], [tl(2),br(2)],'-xr') ;
        plot([tr(1),bl(1)], [tr(2),bl(2)],'-xr') ;
        plot([t(1),b(1)], [t(2),b(2)],'-xr') ;
        plot([l(1),r(1)], [l(2),r(2)],'-xr') ;
        plot([tl12(1),tr12(1)], [tl12(2),tr12(2)],'-xc') ;
        plot([tl12(1),bl12(1)], [tl12(2),bl12(2)],'-xc') ;
        plot([tr12(1),br12(1)], [tr12(2),br12(2)],'-xc') ;
        plot([bl12(1),br12(1)], [bl12(2),br12(2)],'-xc') ;
        plot([tl12(1),br12(1)], [tl12(2),br12(2)],'-xc') ;
        plot([tr12(1),bl12(1)], [tr12(2),bl12(2)],'-xc') ;
        plot([t12(1),b12(1)], [t12(2),b12(2)],'-xc') ;
        plot([l12(1),r12(1)], [l12(2),r12(2)],'-xc') ;
        hold off ;
    end

    txmat = [1 2 3; 1 2 3; 1 2 3] ;
    tymat = txmat' ;

    tdist = zeros(3,3,3,3) ; tang = zeros(3,3,3,3) ;
    for i = 1:3
        for j = 1:3
            for k = 1:3
                for l = 1:3
                    tdist(i,j,k,l) = sqrt( (txmat(k,l) - txmat(i,j)).^2 + (tymat(k,l) - tymat(i,j)).^2 ) ;
                    if ( sqrt( (txmat(k,l) - txmat(i,j)).^2 + (tymat(k,l) - tymat(i,j)).^2 ) == NaN )
                        tdist(i,j,k,l) = 0 ;
                    end
                    tang(i,j,k,l) = atan2( (tymat(k,l) - tymat(i,j)),(txmat(k,l) - txmat(i,j)) ) ;
                    if ( atan2( (tymat(k,l) - tymat(i,j)),(txmat(k,l) - txmat(i,j)) ) == NaN )
                        tang(i,j,k,l) = 0 ;
                    end
                end
            end
        end
    end

    ppxco = [tl12(1) t12(1) tr12(1); l12(1) cent12(1) r12(1); bl12(1) b12(1) br12(1)] ;
    ppyco = [tl12(2) t12(2) tr12(2); l12(2) cent12(2) r12(2); bl12(2) b12(2) br12(2)] ;

    pdist = zeros(3,3,3,3) ; pang = zeros(3,3,3,3) ;
    for i = 1:3
        for j = 1:3
            for k = 1:3
                for l = 1:3
                    pdist(i,j,k,l) = sqrt( (ppxco(k,l) - ppxco(i,j)).^2 + (ppyco(k,l) - ppyco(i,j)).^2 ) ;
                    if ( sqrt( (ppxco(k,l) - ppxco(i,j)).^2 + (ppyco(k,l) - ppyco(i,j)).^2 ) == NaN )
                        pdist(i,j,k,l) = 0 ;
                    end
                    pang(i,j,k,l) = atan2( (ppyco(k,l) - ppyco(i,j)),(ppxco(k,l) - ppxco(i,j)) ) ;
                    if ( atan2( (ppyco(k,l) - ppyco(i,j)),(ppxco(k,l) - ppxco(i,j)) ) == NaN )
                        pang(i,j,k,l) = 0 ;
                    end
                end
            end
        end
    end

    pnormdist = pdist./tdist ;
    redist = reshape(pnormdist,81,1,1,1) ;
    redist = redist(~isnan(redist)) ;
    meandist = mean(redist) ;
    stddist = std(redist) ;

    pnormang = pang - tang ;
    reang = reshape(pnormang,81,1,1,1) ;
    reang = reang(~isnan(reang)) ;
    reang2 = reang ;
    reang2(reang < -1 ) = reang(reang < -1 ) + (2.*pi) ;
    meanang = mean(reang2) ;
    stdang = std(reang2) ;

    gridwidth = 2.*meandist ;
    sp2sp = gridwidth./40 ;

    [spotmat] = ppmakegrid41(cent12,sp2sp,meanang) ;
%    [spotmat] = ppmakegrid41(cent12,sp2sp,0) ;

%    figure, imagesc(ppimi(:,:),[ppmin ppmax]), axis image, hold on ;
%    colormap(gray) ;
    figure, hold on ;
    plot(cent(1),cent(2),'xr') ;    % Mark central point
    plot(cent12(1),cent12(2),'xc') ;    % Mark new central point
    plot([tl(1),tr(1)], [tl(2),tr(2)],'-xr') ;
    plot([tl(1),bl(1)], [tl(2),bl(2)],'-xr') ;
    plot([tr(1),br(1)], [tr(2),br(2)],'-xr') ;
    plot([bl(1),br(1)], [bl(2),br(2)],'-xr') ;
    plot([tl(1),br(1)], [tl(2),br(2)],'-xr') ;
    plot([tr(1),bl(1)], [tr(2),bl(2)],'-xr') ;
    plot([t(1),b(1)], [t(2),b(2)],'-xr') ;
%    plot([l(1),r(1)], [l(2),r(2)],'-xr') ;
    plot([tl12(1),tr12(1)], [tl12(2),tr12(2)],'-xc') ;
    plot([tl12(1),bl12(1)], [tl12(2),bl12(2)],'-xc') ;
    plot([tr12(1),br12(1)], [tr12(2),br12(2)],'-xc') ;
    plot([bl12(1),br12(1)], [bl12(2),br12(2)],'-xc') ;
    plot([tl12(1),br12(1)], [tl12(2),br12(2)],'-xc') ;
    plot([tr12(1),bl12(1)], [tr12(2),bl12(2)],'-xc') ;
    plot([t12(1),b12(1)], [t12(2),b12(2)],'-xc') ;
    plot([l12(1),r12(1)], [l12(2),r12(2)],'-xc') ;
    for i = 1:41
        for j = 1:41
            plot(spotmat(i,j,1),spotmat(i,j,2),'xc') ;
        end
    end
    axis image ;
    hold off ;

%    Refind spot positions

    newspotmat = zeros(41,41,2) ;

    for i = 1:41
        for j = 1:41
            spotpos = round([spotmat(i,j,1) spotmat(i,j,2)]) ;
            [xypos,stdxy] = ppspot(ppim,spotpos,spotwidth,pplim) ;
            newspotmat(i,j,:) = xypos ;
        end
    end

    if 0
        figure, imscale(ppim(:,:),[ppmin ppmax]), axis image, hold on ;
        colormap(gray) ;
        plot(cent(1),cent(2),'xr') ;    % Mark central point
        for i = 1:41
            for j = 1:41
                plot(spotmat(i,j,1),spotmat(i,j,2),'xr') ;
                plot(newspotmat(i,j,1),newspotmat(i,j,2),'xy') ;
            end
        end
        hold off
    end

    ppcent = newspotmat ;
    clear newspotmat ;

%    Calculate distance to nearest neighbour
%    ppnearest matrix is ordered as follows:
%    Horz position
%    Vert position
%    Horz distance to RH neighbour
%    Vert distance to RH neighbour
%    Distance to RH neighbour
%    Angle to RH neighbour
%    Horz distance to lower neighbour
%    Vert distance to lower neighbour
%    Distance to lower neighbour
%    Angle to lower neighbour

    ppnearest = zeros(41,41,10) ;
    ppnearest(:,:,1:2) = ppcent ;

    for i = 1:41
        for j = 1:41
            if j == 41
                ppnearest(i,j,3:6) = 0 ;
            else
                ppnearest(i,j,3:4) = ( ppcent(i,j+1,1:2) - ppcent(i,j,1:2) ) ;
                ppnearest(i,j,5) = sqrt( ppnearest(i,j,3).^2 + ppnearest(i,j,4).^2 ) ;
                ppnearest(i,j,6) = atan( ppnearest(i,j,4)./ppnearest(i,j,3) ) ;
            end
            if i == 41
                ppnearest(i,j,7:10) = 0 ;
            else
                ppnearest(i,j,7:8) = ( ppcent(i+1,j,1:2) - ppcent(i,j,1:2) ) ;
                ppnearest(i,j,9) = sqrt( ppnearest(i,j,7).^2 + ppnearest(i,j,8).^2 ) ;
                ppnearest(i,j,10) = atan( ppnearest(i,j,7)./ppnearest(i,j,8) ) ;
            end
        end
    end    

    tlbr = sqrt( ( ppnearest(41,41,1) - ppnearest(1,1,1) ).^2 + ( ppnearest(41,41,2) - ppnearest(1,1,2) ).^2 ) ;
    trbl = sqrt( ( ppnearest(41,1,1) - ppnearest(1,41,1) ).^2 + ( ppnearest(41,1,2) - ppnearest(1,41,2) ).^2 ) ;
    tltr = sqrt( ( ppnearest(1,41,1) - ppnearest(1,1,1) ).^2 + ( ppnearest(1,41,2) - ppnearest(1,1,2) ).^2 ) ;
    blbr = sqrt( ( ppnearest(41,41,1) - ppnearest(41,1,1) ).^2 + ( ppnearest(41,41,2) - ppnearest(41,1,2) ).^2 ) ;

    ppdistx = mean(reshape(ppnearest(:,1:40,5),(41.*40),1)) ;
    ppangx = mean(reshape(ppnearest(:,1:40,6),(41.*40),1)) ;
    ppdisty = mean(reshape(ppnearest(1:40,:,9),(41.*40),1)) ;
    ppangy = mean(reshape(ppnearest(1:40,:,10),(41.*40),1)) ;
    ppdistxstd = std(reshape(ppnearest(:,1:40,5),(41.*40),1)) ;
    ppangxstd = std(reshape(ppnearest(:,1:40,6),(41.*40),1)) ;
    ppdistystd = std(reshape(ppnearest(1:40,:,9),(41.*40),1)) ;
    ppangystd = std(reshape(ppnearest(1:40,:,10),(41.*40),1)) ;

%    Calculate distance to centre
%    ppcentdist matrix is ordered as follows:
%    Horz position
%    Vert position
%    Horz distance to centre
%    Vert distance to centre
%    Distance to centre
%    Angle to centre
%    Light intensity
%    Normalised light intensity (compared to central spot)

    ppcentdist = zeros(41,41,8) ;
    ppcentdist(:,:,1:2) = ppcent ;

    for i = 1:41
        for j = 1:41
            if ( ( j == 21 ) & ( i == 21 ) )
                ppcentdist(i,j,3:6) = 0 ;
            else
                ppcentdist(i,j,3:4) = ( ppcent(i,j,1:2) - ppcent(21,21,1:2) ) ;
                ppcentdist(i,j,5) = sqrt( ppcentdist(i,j,3).^2 + ppcentdist(i,j,4).^2 ) ;
                ppcentdist(i,j,6) = atan2( ppcentdist(i,j,4),ppcentdist(i,j,3) ) ;
            end
            spotpos = [ppcentdist(i,j,1) ppcentdist(i,j,2)] ;
            [xypos,stdxy,spint] = ppspot(ppim,spotpos,intspotw,pplim) ;
            ppcentdist(i,j,7) = spint ;
        end
    end

    ppcentdist(:,:,8) = ppcentdist(:,:,7)./ppcentdist(21,21,7) ;

    ppcalmat = zeros(41,41,3) ;

    meandists = ppcentdist(:,:,5) ;
    meanangs = ppcentdist(:,:,6) - ppangx ;

    ppcalmat(:,:,2) = ppcentdist(:,:,6) - ppangx ;

    figure, bar3(ppcentdist(:,:,8)) ;

    ppcalmat(:,:,3) = ppcentdist(:,:,7)./ppcentdist(21,21,7) ;
    figure, bar3(ppcalmat(:,:,3)) ;

    figure, hold on ;
    for i = 1:41
        for j = 1:41
            xco = ppcentdist(i,j,5).*cos(ppcentdist(i,j,6)) ;
            yco = ppcentdist(i,j,5).*sin(ppcentdist(i,j,6)) ;
            plot(xco,yco,'xr') ;
            xco2 = ppcent(i,j,1) - ppcent(21,21,1) ;
            yco2 = ppcent(i,j,2) - ppcent(21,21,2) ;
            plot(xco2,yco2,'xc') ;
        end
    end
    hold off

    figure, hold on ;
    for i = 1:41
        for j = 1:41
            xco = meandists(i,j).*cos(meanangs(i,j)) ;
            yco = meandists(i,j).*sin(meanangs(i,j)) ;
            plot(xco,yco,'xr') ;
            xco2 = ppcent(i,j,1) - ppcent(21,21,1) ;
            yco2 = ppcent(i,j,2) - ppcent(21,21,2) ;
            plot(xco2,yco2,'xc') ;
        end
    end
    hold off

    ppspacing = ( ppdistx + ppdisty )./2 ;
    sppix = ppspacing./3 ;   % Number of pixels per mm for spot image
%    load E:\pepperpot\cal_data\rupix.mat ;
%    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 ;
%    midlen = ppspacing.*40 ;

    ppcalmat(:,:,1) = meandists./ppspacing ;

    save('E:\pepperpot\cal_data\ppmk2_caldat.mat','ppcalmat') ;

    break

%    clear ppcalmat ;
%    load E:\pepperpot\cal_data\ppcaldat_good.mat ;
%    load E:\pepperpot\cal_data\pixcal.mat ;
%    rerulen = midlen./mlrulratio ;

%    Recalculate spot positions for each image

    mid = [ppcent(21,21,1) ppcent(21,21,2)] ;
    [ppimreco] = ppspotreco(ppcalmat,mid,rerulen,ppangx,invperscale) ;

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

    return


