% PPCALIB
%
%    PPCALIB - Calibrate pepperpot images

    clear all ; close all ;

    [ppim,stdm,pphist,ppminmax] = pepperplotfun('E:\pepperpot\06_04_06\4ms_10MHz_2ADC_f8_cam-back_50.tif',0,[2 0]) ;

%    ppmin = min(min(ppim)) ; ppmax = max(max(ppim)) ;
    ppmin = ppminmax(1) ; ppmax = ppminmax(2) ;

    cent = [982 1055] ;
    tl = [253 279] ;
    tr = [1764 320] ;
    bl = [197 1792] ;
    br = [1711 1833] ;
    r = [1738 1076] ;
    l = [224 1034] ;
    t = [1007 299] ;
    b = [953 1812] ;

    figure, imagesc(ppim(:,:),[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(cent(1),cent(2),'xr') ;   % Mark central point
%    axis([970 1026 980 1034]) ;
    hold off ;

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

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

    figure, imagesc(ppcut2), axis image, hold on ;
    colormap(gray) ;
    axis tight ;
    hold off ;

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

    [cent2,cent2std] = ppspot(ppim,cent,spotwidth,pplim) ;
    [br2,br2std] = ppspot(ppim,br,spotwidth,pplim) ;
    [bl2,bl2std] = ppspot(ppim,bl,spotwidth,pplim) ;
    [tr2,tr2std] = ppspot(ppim,tr,spotwidth,pplim) ;
    [tl2,tl2std] = ppspot(ppim,tl,spotwidth,pplim) ;
    [b2,b2std] = ppspot(ppim,b,spotwidth,pplim) ;
    [l2,l2std] = ppspot(ppim,l,spotwidth,pplim) ;
    [r2,r2std] = ppspot(ppim,r,spotwidth,pplim) ;
    [t2,t2std] = ppspot(ppim,t,spotwidth,pplim) ;

    figure, imagesc(ppim(:,:),[1 16384]), 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(cent2(1),cent2(2),'xc') ;    % Mark new central point
    plot(tl2(1),tl2(2),'xc') ;      % Mark new upper left point
    plot(tr2(1),tr2(2),'xc') ;      % Mark new upper right point
    plot(bl2(1),bl2(2),'xc') ;      % Mark new lower left point
    plot(br2(1),br2(2),'xc') ;      % Mark new lower right point
    plot(l2(1),l2(2),'xc') ;          % Mark new left point
    plot(t2(1),t2(2),'xc') ;          % Mark new top point
    plot(b2(1),b2(2),'xc') ;          % Mark new bottom point
    plot(r2(1),r2(2),'xc') ;          % Mark new right point
    hold off ;

    figure, imagesc(ppim(:,:),[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(cent(1),cent(2),'xr') ;    % Mark central point
    plot(cent2(1),cent2(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([tl2(1),tr2(1)], [tl2(2),tr2(2)],'-xc') ;
    plot([tl2(1),bl2(1)], [tl2(2),bl2(2)],'-xc') ;
    plot([tr2(1),br2(1)], [tr2(2),br2(2)],'-xc') ;
    plot([bl2(1),br2(1)], [bl2(2),br2(2)],'-xc') ;
    plot([tl2(1),br2(1)], [tl2(2),br2(2)],'-xc') ;
    plot([tr2(1),bl2(1)], [tr2(2),bl2(2)],'-xc') ;
    plot([t2(1),b2(1)], [t2(2),b2(2)],'-xc') ;
    plot([l2(1),r2(1)], [l2(2),r2(2)],'-xc') ;
    hold off ;

    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) = atan( (tymat(k,l) - tymat(i,j))./(txmat(k,l) - txmat(i,j)) ) ;
                    if ( atan( (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 = [tl2(1) t2(1) tr2(1); l2(1) cent2(1) r2(1); bl2(1) b2(1) br2(1)] ;
    ppyco = [tl2(2) t2(2) tr2(2); l2(2) cent2(2) r2(2); bl2(2) b2(2) br2(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) = atan( (ppyco(k,l) - ppyco(i,j))./(ppxco(k,l) - ppxco(i,j)) ) ;
                    if ( atan( (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 < 0 ) = reang(reang < 0 ) + pi ;
    meanang = mean(reang2) ;
    stdang = std(reang2) ;

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

    [spotmat] = ppmakegrid(cent2,sp2sp,meanang) ;

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

%    Refind spot positions

    newspotmat = zeros(25,25,2) ;

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

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

    newspotmat2 = zeros(25,25,2) ;

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

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

    ppcent = newspotmat2 ;
    clear newspotmat newspotmat2 ;

%    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(25,25,10) ;
    ppnearest(:,:,1:2) = ppcent ;

    for i = 1:25
        for j = 1:25
            if j == 25
                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 == 25
                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(25,25,1) - ppnearest(1,1,1) ).^2 + ( ppnearest(25,25,2) - ppnearest(1,1,2) ).^2 ) ;
    trbl = sqrt( ( ppnearest(25,1,1) - ppnearest(1,25,1) ).^2 + ( ppnearest(25,1,2) - ppnearest(1,25,2) ).^2 ) ;
    tltr = sqrt( ( ppnearest(1,25,1) - ppnearest(1,1,1) ).^2 + ( ppnearest(1,25,2) - ppnearest(1,1,2) ).^2 ) ;
    blbr = sqrt( ( ppnearest(25,25,1) - ppnearest(25,1,1) ).^2 + ( ppnearest(25,25,2) - ppnearest(25,1,2) ).^2 ) ;

    ppdistx = mean(reshape(ppnearest(:,1:24,5),(25.*24),1)) ;
    ppangx = mean(reshape(ppnearest(:,1:24,6),(25.*24),1)) ;
    ppdisty = mean(reshape(ppnearest(1:24,:,9),(25.*24),1)) ;
    ppangy = mean(reshape(ppnearest(1:24,:,10),(25.*24),1)) ;
    ppdistxstd = std(reshape(ppnearest(:,1:24,5),(25.*24),1)) ;
    ppangxstd = std(reshape(ppnearest(:,1:24,6),(25.*24),1)) ;
    ppdistystd = std(reshape(ppnearest(1:24,:,9),(25.*24),1)) ;
    ppangystd = std(reshape(ppnearest(1:24,:,10),(25.*24),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(25,25,8) ;
    ppcentdist(:,:,1:2) = ppcent ;

    for i = 1:25
        for j = 1:25
            if ( ( j == 13 ) & ( i == 13 ) )
                ppcentdist(i,j,3:6) = 0 ;
            else
                ppcentdist(i,j,3:4) = ( ppcent(i,j,1:2) - ppcent(13,13,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,pplimlow) ;
            ppcentdist(i,j,7) = spint ;
        end
    end

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

%    Compare to another image

    [ppim2,stdm2,pphist2,ppminmax2] = pepperplotfun('E:\pepperpot\06_04_06\4ms_10MHz_2ADC_f8_cam-in_50.tif',0,[2 0]) ;
    centpp2 = [1053 1106] ;
    ppl = [298 1086] ;
    ppr = [1810 1127] ;

    figure, imagesc(ppim2,[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(centpp2(1),centpp2(2),'xr') ;   % Mark central point
    plot(ppl(1),ppl(2),'xr') ;
    plot(ppr(1),ppr(2),'xr') ;
    hold off ;

    [ppspotmat] = ppmakegrid(centpp2,sp2sp,meanang) ;

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

%    Refind spot positions

    newspotmat = zeros(25,25,2) ;

    for i = 1:25
        for j = 1:25
            spotpos = round([ppspotmat(i,j,1) ppspotmat(i,j,2)]) ;
            [xypos,stdxy] = ppspot(ppim2,spotpos,spotwidth,pplim) ;
            newspotmat(i,j,:) = xypos ;
        end
    end

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

    newspotmat2 = zeros(25,25,2) ;

    for i = 1:25
        for j = 1:25
            spotpos = [newspotmat(i,j,1) newspotmat(i,j,2)] ;
            [xypos,stdxy] = ppspot(ppim2,spotpos,spotwidth,pplim) ;
            newspotmat2(i,j,:) = xypos ;
        end
    end

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

    ppcent2 = newspotmat2 ;
    clear newspotmat newspotmat2 ;

%    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

    ppnearest2 = zeros(25,25,10) ;
    ppnearest2(:,:,1:2) = ppcent2 ;

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

    tlbr2 = sqrt( ( ppnearest2(25,25,1) - ppnearest2(1,1,1) ).^2 + ( ppnearest2(25,25,2) - ppnearest2(1,1,2) ).^2 ) ;
    trbl2 = sqrt( ( ppnearest2(25,1,1) - ppnearest2(1,25,1) ).^2 + ( ppnearest2(25,1,2) - ppnearest2(1,25,2) ).^2 ) ;
    tltr2 = sqrt( ( ppnearest2(1,25,1) - ppnearest2(1,1,1) ).^2 + ( ppnearest2(1,25,2) - ppnearest2(1,1,2) ).^2 ) ;
    blbr2 = sqrt( ( ppnearest2(25,25,1) - ppnearest2(25,1,1) ).^2 + ( ppnearest2(25,25,2) - ppnearest2(25,1,2) ).^2 ) ;

    ppdistx2 = mean(reshape(ppnearest2(:,1:24,5),(25.*24),1)) ;
    ppangx2 = mean(reshape(ppnearest2(:,1:24,6),(25.*24),1)) ;
    ppdisty2 = mean(reshape(ppnearest2(1:24,:,9),(25.*24),1)) ;
    ppangy2 = mean(reshape(ppnearest2(1:24,:,10),(25.*24),1)) ;
    ppdistx2std = std(reshape(ppnearest2(:,1:24,5),(25.*24),1)) ;
    ppangx2std = std(reshape(ppnearest2(:,1:24,6),(25.*24),1)) ;
    ppdisty2std = std(reshape(ppnearest2(1:24,:,9),(25.*24),1)) ;
    ppangy2std = std(reshape(ppnearest2(1:24,:,10),(25.*24),1)) ;

%    Calculate distance to centre
%    ppcentdist2 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)

    ppcentdist2 = zeros(25,25,8) ;
    ppcentdist2(:,:,1:2) = ppcent2 ;

    for i = 1:25
        for j = 1:25
            if ( ( j == 13 ) & ( i == 13 ) )
                ppcentdist2(i,j,3:6) = 0 ;
            else
                ppcentdist2(i,j,3:4) = ( ppcent2(i,j,1:2) - ppcent2(13,13,1:2) ) ;
                ppcentdist2(i,j,5) = sqrt( ppcentdist2(i,j,3).^2 + ppcentdist2(i,j,4).^2 ) ;
                ppcentdist2(i,j,6) = atan2( ppcentdist2(i,j,4),ppcentdist2(i,j,3) ) ;
            end
            spotpos = [ppcentdist2(i,j,1) ppcentdist2(i,j,2)] ;
            [xypos,stdxy,spint] = ppspot(ppim2,spotpos,intspotw,pplimlow) ;
            ppcentdist2(i,j,7) = spint ;
        end
    end

    ppcentdist2(:,:,8) = ppcentdist2(:,:,7)./ppcentdist2(13,13,7) ;

%    Compare to third image

    [ppim3,stdm3,pphist3,ppminmax3] = pepperplotfun('E:\pepperpot\06_04_06\4ms_10MHz_2ADC_f8_cam-mid_50.tif',0,[2 0]) ;
    centpp3 = [1033 1081] ;
    ppl3 = [278 1060] ;
    ppr3 = [1790 1102] ;

    figure, imagesc(ppim3,[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(centpp3(1),centpp3(2),'xr') ;   % Mark central point
    plot(ppl3(1),ppl3(2),'xr') ;
    plot(ppr3(1),ppr3(2),'xr') ;
    hold off ;

    [ppspotmat3] = ppmakegrid(centpp3,sp2sp,meanang) ;

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

%    Refind spot positions

    newspotmat = zeros(25,25,2) ;

    for i = 1:25
        for j = 1:25
            spotpos = round([ppspotmat3(i,j,1) ppspotmat3(i,j,2)]) ;
            [xypos,stdxy] = ppspot(ppim3,spotpos,spotwidth,pplim) ;
            newspotmat(i,j,:) = xypos ;
        end
    end

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

    newspotmat2 = zeros(25,25,2) ;

    for i = 1:25
        for j = 1:25
            spotpos = [newspotmat(i,j,1) newspotmat(i,j,2)] ;
            [xypos,stdxy] = ppspot(ppim3,spotpos,spotwidth,pplim) ;
            newspotmat2(i,j,:) = xypos ;
        end
    end

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

    ppcent3 = newspotmat2 ;
    clear newspotmat newspotmat2 ;

%    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

    ppnearest3 = zeros(25,25,10) ;
    ppnearest3(:,:,1:2) = ppcent3 ;

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

    tlbr3 = sqrt( ( ppnearest3(25,25,1) - ppnearest3(1,1,1) ).^2 + ( ppnearest3(25,25,2) - ppnearest3(1,1,2) ).^2 ) ;
    trbl3 = sqrt( ( ppnearest3(25,1,1) - ppnearest3(1,25,1) ).^2 + ( ppnearest3(25,1,2) - ppnearest3(1,25,2) ).^2 ) ;
    tltr3 = sqrt( ( ppnearest3(1,25,1) - ppnearest3(1,1,1) ).^2 + ( ppnearest3(1,25,2) - ppnearest3(1,1,2) ).^2 ) ;
    blbr3 = sqrt( ( ppnearest3(25,25,1) - ppnearest3(25,1,1) ).^2 + ( ppnearest3(25,25,2) - ppnearest3(25,1,2) ).^2 ) ;

    ppdistx3 = mean(reshape(ppnearest3(:,1:24,5),(25.*24),1)) ;
    ppangx3 = mean(reshape(ppnearest3(:,1:24,6),(25.*24),1)) ;
    ppdisty3 = mean(reshape(ppnearest3(1:24,:,9),(25.*24),1)) ;
    ppangy3 = mean(reshape(ppnearest3(1:24,:,10),(25.*24),1)) ;
    ppdistx3std = std(reshape(ppnearest3(:,1:24,5),(25.*24),1)) ;
    ppangx3std = std(reshape(ppnearest3(:,1:24,6),(25.*24),1)) ;
    ppdisty3std = std(reshape(ppnearest3(1:24,:,9),(25.*24),1)) ;
    ppangy3std = std(reshape(ppnearest3(1:24,:,10),(25.*24),1)) ;

%    Calculate distance to centre
%    ppcentdist3 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)

    ppcentdist3 = zeros(25,25,8) ;
    ppcentdist3(:,:,1:2) = ppcent3 ;

    for i = 1:25
        for j = 1:25
            if ( ( j == 13 ) & ( i == 13 ) )
                ppcentdist3(i,j,3:6) = 0 ;
            else
                ppcentdist3(i,j,3:4) = ( ppcent3(i,j,1:2) - ppcent3(13,13,1:2) ) ;
                ppcentdist3(i,j,5) = sqrt( ppcentdist3(i,j,3).^2 + ppcentdist3(i,j,4).^2 ) ;
                ppcentdist3(i,j,6) = atan2( ppcentdist3(i,j,4),ppcentdist3(i,j,3) ) ;
            end
            spotpos = [ppcentdist3(i,j,1) ppcentdist3(i,j,2)] ;
            [xypos,stdxy,spint] = ppspot(ppim3,spotpos,intspotw,pplimlow) ;
            ppcentdist3(i,j,7) = spint ;
        end
    end

    ppcentdist3(:,:,8) = ppcentdist3(:,:,7)./ppcentdist3(13,13,7) ;

    ppcut3 = ( ppim3 > pplim ) ;
    ppimcut3 = ( ppim3 - pplim ).*ppcut3 ;
    ppcut3a = ( ppim3 > pplim2 ) ;
    ppimcut3a = ( ppim3 - pplim2 ).*ppcut3a ;
    ppcut3b = ( ppim3 > pplimlow ) ;
    ppimcut3b= ( ppim3 - pplimlow ).*ppcut3b ;

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

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

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

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

    [ppspotmat3a] = round(ppmakegrid(squeeze(ppcent3(13,13,:)),ppdistx3,ppangx3,ppdisty3,ppangy3)) ;

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

    figure, imagesc(ppim3(:,:),[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(centpp3(1),centpp3(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            xco = ppcentdist3(i,j,3) ;
            yco = ppcentdist3(i,j,4) ;
            plot((xco+ppcentdist3(13,13,1)),(yco+ppcentdist3(13,13,2)),'xc') ;
            xco2 = ppcentdist3(i,j,5).*cos(ppcentdist3(i,j,6)) ;
            yco2 = ppcentdist3(i,j,5).*sin(ppcentdist3(i,j,6)) ;
            plot((xco2+ppcentdist3(13,13,1)),(yco2+ppcentdist3(13,13,2)),'xm') ;
        end
    end
    hold off

    figure, imagesc(ppim3(:,:),[1 16384]), axis image, hold on ;
    colormap(gray) ;
    plot(centpp3(1),centpp3(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            xco = ppcentdist(i,j,5).*cos(ppcentdist(i,j,6)-ppangx) ;
            yco = ppcentdist(i,j,5).*sin(ppcentdist(i,j,6)-ppangx) ;
            plot((xco+ppcentdist3(13,13,1)),(yco+ppcentdist3(13,13,2)),'xc') ;
            xco = ppcentdist2(i,j,5).*cos(ppcentdist2(i,j,6)-ppangx2) ;
            yco = ppcentdist2(i,j,5).*sin(ppcentdist2(i,j,6)-ppangx2) ;
            plot((xco+ppcentdist3(13,13,1)),(yco+ppcentdist3(13,13,2)),'xm') ;
            xco = ppcentdist3(i,j,5).*cos(ppcentdist3(i,j,6)-ppangx3) ;
            yco = ppcentdist3(i,j,5).*sin(ppcentdist3(i,j,6)-ppangx3) ;
            plot((xco+ppcentdist3(13,13,1)),(yco+ppcentdist3(13,13,2)),'xy') ;
        end
    end
    hold off

    ppcalmat = zeros(25,25,3) ;

    alldists = zeros(25,25,3) ;
    alldists(:,:,1) = ppcentdist(:,:,5) ;
    alldists(:,:,2) = ppcentdist2(:,:,5) ;    
    alldists(:,:,3) = ppcentdist3(:,:,5) ;
    meandists = mean(alldists,3) ;

    allangs = zeros(25,25,3) ;
    allangs(:,:,1) = ppcentdist(:,:,6) - ppangx ;
    allangs(:,:,2) = ppcentdist2(:,:,6) - ppangx2 ;
    allangs(:,:,3) = ppcentdist3(:,:,6) - ppangx3 ;
    meanangs = mean(allangs,3) ;

    ppcalmat(:,:,2) = meanangs ;

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

    alllightint = zeros(25,25,3) ;
    alllightint(:,:,1) = ppcentdist(:,:,7) ;
    alllightint(:,:,2) = ppcentdist2(:,:,7) ;
    alllightint(:,:,3) = ppcentdist3(:,:,7) ;
    meanlightint = mean(alllightint,3) ;

    ppcalmat(:,:,3) = meanlightint./meanlightint(13,13) ;
    figure, bar3(ppcalmat(:,:,3)) ;

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            xco = meandists(i,j).*cos(meanangs(i,j)) ;
            yco = meandists(i,j).*sin(meanangs(i,j)) ;
            plot(xco,yco,'xr') ;
        end
    end
    hold off

    ppspacing = ( ppdistx + ppdistx2 + ppdistx3 + ppdisty + ppdisty2 + ppdisty3 )./6 ;
%    sppix = ( ppspacing .* 24 )./72 ;   % 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.*24 ;

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

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

    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(13,13,1) ppcent(13,13,2)] ;
    [ppimreco] = ppspotreco(ppcalmat,mid,rerulen,ppangx,invperscale) ;

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

    mid2 = [ppcent2(13,13,1) ppcent2(13,13,2)] ;
    [ppim2reco] = ppspotreco(ppcalmat,mid2,rerulen,ppangx2,invperscale) ;

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

    mid3 = [ppcent3(13,13,1) ppcent3(13,13,2)] ;
    [ppim3reco] = ppspotreco(ppcalmat,mid3,rerulen,ppangx3,invperscale) ;

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

    return


