clear all
close all

%% Here we set the parameters for the ploynomial fit and the plotting part

Font_Size = 13;  % set the size of the typo
rs = 600;  % Set the height range of the height map image
n = 1; % set the degre of the fiting polynomial
ps = 87.89;  %This is equal to the scan size divided by the number of pixels, it's in nm
m = 4; % scaling parameter for the plots

%% Load the cropped images with Gwyddion

Z1 = load('A4_35x35um_pre_test')
Z2 = load('A4_35x35um_post_test')

%% Here we put the height maps in nm

Z_pre_test_norot = Z1.x90A4*10^9
Z_post_test_norot = Z2.A435x35croppedposttest*10^9

Z_pre_test_norot_reshaped = reshape(Z_pre_test_norot, length(Z_pre_test_norot(:,1))*length(Z_pre_test_norot(1,:)),1);
rms_h_rough5s_3bar = std(Z_pre_test_norot_reshaped)

%% Here we rotate one image to the other by 10 degres, tis value will differ depending on the data

Z_pre_test = imrotate(Z_pre_test_norot, 0)
Z_post_test = imrotate(Z_post_test_norot, -1) 

%% Here we shrink one image to the other form the center of each image, such that they have the same size

center1 = [ceil(numel(Z_pre_test(1,:))/2),ceil(numel(Z_pre_test(:,1))/2)]; %select a centering region for the pre-test surface to overlap the two data afterhand
center2 = [ceil(numel(Z_post_test(1,:))/2),ceil(numel(Z_post_test(:,1))/2)]; %select a centering region for the post-test surface to overlap the two data afterhand
range = 170; % select how large you want the region to be: + and - the range in both axis

Z_pre_test_centered_ = Z_pre_test(center1(1)-range:center1(1)+range,center1(2)-range:center1(2)+range); %% define the new image size
Z_pre_test_centered = Z_pre_test_centered_-mean(mean(Z_pre_test_centered_))+rs; %data manipulation before plotting it


Z_post_test_centered_ = Z_post_test(center2(1)-range:center2(1)+range,center2(2)-range:center2(2)+range); %% define the new image size
Z_post_test_centered = Z_post_test_centered_-mean(mean(Z_post_test_centered_))+rs; %data manipulation before plotting it

% Z_pre_test_fit = Z_pre_test_centered;
% Z_post_test_fit = Z_post_test_centered;

%% Here we fit the new height map with a polynomial of degree 1
% 
x = (1:length(Z_pre_test_centered))*ps;
y = (1:length(Z_post_test_centered))*ps;
 
for i=1:length(Z_pre_test_centered(:,1))
    p = polyfit(x,Z_pre_test_centered(i,:),n);
    Z_pre_test_fit(i,:) = Z_pre_test_centered(i,:) - polyval(p,x);
end

 
for i=1:length(Z_post_test_centered(:,1))
    p = polyfit(y,Z_post_test_centered(i,:),n);
    Z_post_test_fit(i,:) = Z_post_test_centered(i,:) - polyval(p,y);
end

%% We plot here the fitted and centered images before and after

figure
image((Z_pre_test_fit-mean(mean(Z_pre_test_fit))+rs)/(2*rs)*64); 
axis off
set(gca,'plotboxaspectratio',[1 1 1])
movegui(gca,'northwest')
colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14)
colormap jet

figure
image((Z_post_test_fit-mean(mean(Z_post_test_fit))+rs)/(2*rs)*64); 
axis off
set(gca,'plotboxaspectratio',[1 1 1])
movegui(gca,'northwest')
colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14)
colormap jet

%% Now we apply the plot_difference functionto find the optimal vertical and horizontal shiftfrom one map to the other

min_overlap = plot_difference(Z_pre_test_fit,Z_post_test_fit,20);

[opt_shift_value, shift_index] = min(min_overlap(:)); % We define here the optimal shit value
[I,J] = ind2sub([size(min_overlap,1) size(min_overlap,2)], shift_index);

opti_index_shift = [I,J];

%% We redefine here a new center and re-crop the pre-height map to have it optimally align with the post-image. We also re-apply the polynomial fit in case of the the alignment procedure influenced the shape of the data

center3 = [center1(1)+J-39, center1(2)+I-35]
range = 125;
center4 = [ceil(numel(Z_post_test(1,:))/2),ceil(numel(Z_post_test(:,1))/2)];
Z_pre_test_opti_ = Z_pre_test(center3(1)-range:center3(1)+range,center3(2)-range:center3(2)+range); %% define the new image size
Z_pre_test_opti = Z_pre_test_opti_-mean(mean(Z_pre_test_opti_))+rs; %data manipulation before plotting it

Z_post_test_opti_ = Z_post_test(center4(1)-range:center4(1)+range,center4(2)-range:center4(2)+range); %% define the new image size
Z_post_test_opti = Z_post_test_opti_-mean(mean(Z_post_test_opti_))+rs; %data manipulation before plotting it

% t = (1:length(Z_pre_test_opti))*ps;
%  
% for i=1:length(Z_pre_test_opti(:,1))
%     p = polyfit(t,Z_pre_test_opti(i,:),n);
%     Z_pre_test_opti_fit(i,:) = Z_pre_test_opti(i,:) - polyval(p,t);
% end

Z_pre_test_opti_fit = Z_pre_test_opti;

%% We plot here the new aligned pre test image

figure
%image(Z)
image((Z_pre_test_opti_fit-mean(mean(Z_pre_test_opti_fit))+rs)/(2*rs)*64); 
% image(Z(261:580,261:580)/(2*rs)*64);
axis off
set(gca,'plotboxaspectratio',[1 1 1])
movegui(gca,'northwest')
colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14)
colormap jet


figure
%image(Z)
image((Z_post_test_opti-mean(mean(Z_post_test_opti))+rs)/(2*rs)*64); 
% image(Z(261:580,261:580)/(2*rs)*64);
axis off
set(gca,'plotboxaspectratio',[1 1 1])
movegui(gca,'northwest')
colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14)
colormap jet

%% We define the difference image here, save the corresponding height maps and plot it

opti_height_shift = optimizing_height_diff(Z_pre_test_opti_fit, Z_post_test_fit, 15)

Z_pre_test_opti_fit = Z_pre_test_opti_fit + opti_height_shift

Z_diff = Z_post_test_opti - Z_pre_test_opti_fit

figure
image((Z_diff-mean(mean(Z_diff))+rs)/(2*rs)*64)
axis off
set(gca,'plotboxaspectratio',[1 1 1])
movegui(gca,'northeast')
colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14) 
colormap jet

Z_diff_abs =  abs(Z_diff)

opti_diff = mean(mean(Z_diff))

Z_diff = Z_diff-mean(mean(Z_diff));
a = Z_diff(:)
% a = reshape(I_conv,[1,npixel*npixel]);
% otsu_th = f*graythresh(im_fluor/max(max(im_fluor)))*max(max(im_fluor));
binwidth = 1;
bins = -rs:binwidth:rs;
bincounts_a = histc(a,bins);
y = bincounts_a;
x = bins;
x(y==0) = [];
y(y==0) = [];

y = transpose(y).*x;
y = y*ps^2;
y = transpose(y);

FigHandle = figure('Position', [500,300,480,480]); % the last two numbers define the size of the figure; good to keep it small. If the size is small and it is still readable it will also be good in a paper/presentation
plot(x,abs(y),'LineWidth', 1.5, 'Color', [0 0 0], 'DisplayName','Height distribution difference image') % specify linewidth and color and name in legend
axis([-300 300 -40 180000000])  % specify how to crop the graph; good to make sure there is a bit of space around the data; if you want to include insets you can edit this so they will fit
set(gca, 'linewidth',1.5,'fontsize',13);  %appropriate line width and font size for the size of the graph
set(gca, 'plotboxaspectratio', [1 1 1]);  % Show a box around the graph
xh = xlabel('Heigth \it{H} (nm)');  %label x-axis
yh = ylabel('Volme of meterial per height (H*ps^2*numb of pix) \it{V} (nm^3)');  %label y-axis
set([xh,yh],'fontsize',14);  %Font size of the axis labels
movegui(gca, 'northeast');  %move graph into right top corner of the screen
% legend('Location','nm');  %produce legend



center_4 = [64,111];
range = 20; % select how large you want the region to be: + and - the range in both axis

Z_diff_flat_part = Z_diff(center_4(1)-range:center_4(1)+range,center_4(2)-range:center_4(2)+range); %% define the new image size
% Z_diff_flat_part_normalized = Z_diff_flat_part-mean(mean(Z_diff_flat_part))+rs; %data manipulation before plotting it
Z_diff_flat_part_normalized = Z_diff_flat_part-mean(mean(Z_diff_flat_part)); %data manipulation before plotting it

save('5srough_3bar_iso_unworn_difference.mat','Z_diff_flat_part_normalized')

Z_diff_flat_part_normalized_reshaped = reshape(Z_diff_flat_part_normalized, length(Z_diff_flat_part_normalized(:,1))*length(Z_diff_flat_part_normalized(1,:)),1);
rms_h_flat_part = mean(std(Z_diff_flat_part_normalized))

% Z_pre_test_opti_fit = Z_pre_test_opti_fit(3:258,3:258);
% Z_post_test_fit = Z_post_test_fit(3:258,3:258);
% Z_diff = Z_diff(3:258,3:258);
% % 
% save('5srough_3bar_pre_post_test_difference.mat','Z_diff')
% % xlswrite('5srough_3bar_pre_test_opti_cropped.xlsx', [Z_pre_test_opti_fit*10^-9])
% dlmwrite('5srough_3bar_pre_test_opti_cropped.txt', Z_pre_test_opti_fit*10^-9)
% % xlswrite('5srough_3bar_post_test_opti_cropped.xlsx', [Z_post_test_fit*10^-9])
% dlmwrite('5srough_3bar_post_test_opti_cropped.txt', Z_post_test_opti*10^-9)
% save('5srough_3bar_pre_test_opti_cropped.mat', 'Z_pre_test_opti_fit')
% save('5srough_3bar_post_test_opti_cropped.mat', 'Z_post_test_opti')

rs = 600

figure
image((Z_diff-mean(mean(Z_diff))+rs)/(2*rs)*64)
axis off
set(gca,'plotboxaspectratio',[1 1 1])
movegui(gca,'northeast')
colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14) 
colormap jet

%% Crop maps for individual asperities isolation and plot the cross section of those asperities

% Asp_2_cropped_before = Z_pre_test_opti_fit(205:245,63:103);
% Asp_2_cropped_after = Z_post_test_fit(205:245,63:103);
% 
% figure
% image((Asp_2_cropped_before-mean(mean(Asp_2_cropped_before))+rs)/(2*rs)*64)
% axis off
% set(gca,'plotboxaspectratio',[1 1 1])
% movegui(gca,'northeast')
% colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14) 
% colormap jet
% 
% figure
% image((Asp_2_cropped_after-mean(mean(Asp_2_cropped_after))+rs)/(2*rs)*64)
% axis off
% set(gca,'plotboxaspectratio',[1 1 1])
% movegui(gca,'northeast')
% colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14) 
% colormap jet
% 
% cross_section_before = Asp_2_cropped_before(1:36,29)+21;
% cross_section_after = Asp_2_cropped_after(5:40,29)+41;
% 
% cross_section_length = 1:87.9:3164.4;
% 
% % cross_section_before = Asp_2_cropped_before(23,:)+121;
% % cross_section_after = Asp_2_cropped_after(23,:)+158;
% 
% figure()
% plot(cross_section_length, cross_section_before, 'LineWidth', 1.5, 'Color', [0 0 0], 'DisplayName','Initial height binned');
% hold on;
% plot(cross_section_length, cross_section_after, 'LineWidth', 1.5, 'Color', [1 0 0], 'DisplayName','Initial height binned');
% hold on;
% axis([-1 3200 -250 700])
% set(gca, 'linewidth',1.5,'fontsize',13);  %appropriate line width and font size for the size of the graph
% set(gca, 'plotboxaspectratio', [1 1 1]);  % Show a box around the graph
% xh = xlabel('Cross-section length ');  %label x-axis
% yh = ylabel('Asperity height \it{H} (nm)');  %label y-axis
% set([xh,yh],'fontsize',14);  %Font size of the axis labels
% movegui(gca, 'northeast');  %move graph into right top corner of the screen
% % legend('Track 400', 'Track 700', 'Track 900', 'Track 1000', 'Track 1100');
% 
% center_4 = [226,210];
% range = 30; % select how large you want the region to be: + and - the range in both axis
% 
% Z_diff_flat_part = Z_diff(center_4(1)-range:center_4(1)+range,center_4(2)-range:center_4(2)+range); %% define the new image size
% Z_diff_flat_part_normalized = Z_diff_flat_part-mean(mean(Z_diff_flat_part))+rs; %data manipulation before plotting it
% 
% Z_diff_flat_part_normalized_reshaped = reshape(Z_diff_flat_part_normalized, length(Z_diff_flat_part_normalized(:,1))*length(Z_diff_flat_part_normalized(1,:)),1);
% rms_h_flat_part = mean(std(Z_diff_flat_part_normalized))
% 
% figure
% image((Z_diff_flat_part_normalized-mean(mean(Z_diff_flat_part_normalized))+rs)/(2*rs)*64); 
% axis off
% set(gca,'plotboxaspectratio',[1 1 1])
% movegui(gca,'northwest')
% colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14)
% colormap jet

% %% Here we add to the analysis the plane tilt 
% rs = 900;
% 
% out_tilt = optimizing_tilt(Z_pre_test_opti, Z_post_test_centered, -0.1); % We use the tilting function
% [opt_tilt_value, tilt_index] = min(min(out_tilt(:))) ; % We set here the optimal tilt value
% 
% %% We can also do the tilt manually
% 
% a = -1.5; %positive moves center down
% b = -1.0; %negative moves center to left (after)
% 
% 
% for i=1:length(Z_post_test_centered) % this is the loop to correct a tilt, not sure it's that of a need as the data are already flatten
%    for j = 1:length(Z_post_test_centered)
%        Z_post_test_centered_tilt(i,j) = Z_post_test_centered(i,j)+a*i+b*j; % the loop is based on two coefficient that shift up or down each pixels, according to their number.
%    end
% end
% 
% %% We define the new difference height map here and plot it again
% 
% Z_diff2 = Z_post_test_centered_tilt - Z_pre_test_opti
% 
% opti_diff2 = mean(mean(Z_diff2))
% 
% 
% rs = 1000;
% 
% figure
% image((Z_diff2-mean(mean(Z_diff2))+rs)/(2*rs)*64);
% axis off
% set(gca,'plotboxaspectratio',[1 1 1])
% movegui(gca,'northeast')
% colorbar('YTick',(0:round(2*rs/m):2*rs)/(2*rs)*64,'YTickLabel',(-rs:round(2*rs/m):rs),'FontSize',14) 
% colormap jet

%% The functions used in the scipt are here

function out = plot_difference(pre_image,post_image,c) %Build the function that takes as an input Im_fluor: a serie of matrices stacked on top of each other, c is the range over which we try to overlap

tmp = pre_image(:,:);  %select the pre_image matrix
img = post_image(:,:);  %select the post_image matrix

pos1 = [ceil(numel(tmp(1,:))/2),ceil(numel(tmp(:,1))/2)]; % find automatically the x and y coordinates of the common starting point (center because cropped beforehand) for tmp. Can also enter the coordinates manually
pos2 = [ceil(numel(img(1,:))/2),ceil(numel(img(:,1))/2)]; % do the same but for the other image img

pos_1 = pos1;
pos_2 = pos2;
pos = pos_1;

e = round(pos(1)); % rounded x axis pos coordinate
d = round(pos(2)); % rounded y axis pos coordinate

diff = NaN(2*c+1,2*c+1); % define the difference between the two images

for i = -c+d:c+d  % double loop that will successively shift one image to the other and caculate the difference between the two. If stops when it find the position where the difference is minimal.
    for j = -c+e:c+e
        if i>=0&&j>=0
            diff(i+c+1-d,j+c+1-e) = mean(mean(abs(img(i+1:end,j+1:end)-tmp(1:end-i,1:end-j))));
        elseif i<0&&j>=0
            diff(i+c+1-d,j+c+1-e) = mean(mean(abs(img(1:end+i,j+1:end)-tmp(abs(i)+1:end,1:end-j))));
        elseif i>=0&&j<0
            diff(i+c+1-d,j+c+1-e) = mean(mean(abs(img(i+1:end,1:end+j)-tmp(1:end-i,abs(j)+1:end))));
        elseif i<0&&j<0
            diff(i+c+1-d,j+c+1-e) = mean(mean(abs(img(1:end+i,1:end+j)-tmp(abs(i)+1:end,abs(j)+1:end))));
        end

    end
end

out = diff;

end


function out1 = optimizing_tilt(pre_image,post_image,iteration)

a = -0.1:iteration:-2.0
b = -0.1:iteration:-2.0

ita = 1;
itb = 1;


tilt = 1:400;

for h = a
    for g = b
        
        for i=1:length(post_image) % this is the loop to correct a tilt, not sure it's that of a need as the data are already flatten
           for j = 1:length(post_image)
               post_image_tilt(i,j) = post_image(i,j)+h*i+g*j; % the loop is based on two coefficient that shift up or down each pixels, according to their number.
           
           end
        end
        
        difference_tilt = mean(mean(post_image_tilt - pre_image));
        disp(difference_tilt)
        tilt(itb) = difference_tilt;
         
        itb = itb +1
         
    end
    
end

out1 = reshape((tilt),[20,20]);

end
 
function out4 = optimizing_height_diff(pre_image, post_image, shiftval);

rs = 800;
m =4;
c = shiftval;
a_pre = pre_image(:);
a_post = post_image(:);


binwidth = 2;
bins = -rs:binwidth:rs;
bincounts_pre = histc(a_pre,bins);
bincounts_post = histc(a_post,bins);
y_pre = bincounts_pre;
y_post = bincounts_post;
x_pre = bins;
x_post = bins;

diff = NaN(2*c,1); % define the difference between the two distribution

    for i = -c:1:c  % double loop that will sccessively shift one image to the other and caculate the difference between the two. If stops when it find the position where the difference is minimal.
    
            if i<0
                ilow_pre = y_pre;
                ilow_post = y_post;
                ilow_pre = ilow_pre(abs(i)+1:end);
                ilow_post = ilow_post(1:end-abs(i));
                diff(i+c+1) = mean(abs(ilow_post-ilow_pre));
            elseif i==0
                lol_pre = y_pre;
                lol_post = y_post;
                diff(i+c+1) = mean(abs(lol_post-lol_pre));
            elseif i>0
                ihigh_pre = y_pre;
                ihigh_post = y_post;
                ihigh_pre = ihigh_pre(1:end-abs(i));
                ihigh_post = ihigh_post(abs(i)+1:end);
                diff(i+c+1) = mean(abs(ihigh_post-ihigh_pre));

            end
    end

[opt_shift_value, shift_index] = min(diff(:));
optimal_displacement = shift_index-c-1;


out4 = opt_shift_value

end
