clear all
close all


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

Font_Size = 13;
rs = 50;
n = 2;
ps = (87.89)/2;  %This is equal to the scan size divided by the number of pixels, it's in nm
m = 4;

%% Load the cropped images with Gwyddion

Z1 = load('90um_area4_cropped_1212lines_before.mat')
Z2 = load('90um_cropped_1212lines_after.mat')




%% Here we put the height maps in nm

Z_pre_test = Z1.x90umarea4cropped1209linesbefore*10^9
% Z_post_test = Z2.x90umcropped1209linesafter*10^9
Z_post_test_norot = Z2.x90umcropped1209linesafter*10^9

figure
%image(Z)
image((Z_pre_test-mean(mean(Z_pre_test))+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_norot-mean(mean(Z_post_test_norot))+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

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

Z_post_test = imrotate(Z_post_test_norot, -3)


%% 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 = 500; % 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

%% 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

Z_pre_test_opti_reshaped = reshape(Z_pre_test_fit, length(Z_pre_test_fit(:,1))*length(Z_pre_test_fit(1,:)),1);
rms_h_pristine2 = std(Z_pre_test_opti_reshaped)

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

rs = 50;

min_overlap = plot_difference(Z_pre_test_fit,Z_post_test_fit,20);

[opt_shift_value, shift_index] = min(min_overlap(:));
[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+15, center1(2)+I+18];
range = 256;

Z_pre_test_opti_ = Z_pre_test_fit(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_centered_ = Z_post_test_fit(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

%% We plot here the new aligned pre test image

figure
%image(Z)
image((Z_pre_test_opti-mean(mean(Z_pre_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

figure
%image(Z)
image((Z_post_test_centered-mean(mean(Z_post_test_centered))+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
 

%% Here we add to the analysis the plane tilt 

% out_tilt = optimizing_tilt(Z_pre_test_opti, Z_post_test_centered, 0.1)
% [opt_tilt_value, tilt_index] = min(out_tilt(:))

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

opti_height_shift = optimizing_height_diff(Z_pre_test_opti, Z_post_test_centered, 15)

Z_pre_test_opti = Z_pre_test_opti + opti_height_shift

Z_diff = (Z_post_test_centered) - (Z_pre_test_opti)

Z_diff_abs =  abs(Z_diff)

opti_diff = mean(mean(Z_diff))

Z_pre_test_opti = Z_pre_test_opti(1:512,1:512);
Z_post_test_centered = Z_post_test_centered(1:512,1:512);
Z_diff = Z_diff(1:512,1:512);

center_4 = [182,53];
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))

% 
% Z_pre_test_opti_reshaped = reshape(Z_pre_test_opti, length(Z_pre_test_opti(:,1))*length(Z_pre_test_opti(1,:)),1);
% rms_h_pristine2 = std(Z_pre_test_opti_reshaped)

save('pristine2_pre_post_test_difference.mat','Z_diff')
% xlswrite('Pristine2_pre_test_opti_cropped.xlsx', [Z_pre_test_opti*10^-9])
% xlswrite('Pristine2_post_test_opti_cropped.xlsx', [Z_post_test_centered*10^-9])
save('pristine2_pre_test_opti_crop.mat','Z_pre_test_opti')
save('pristine2_post_test_opti_crop.mat','Z_post_test_centered')
dlmwrite('pristine2_pre_test_opti_crop.txt', ((Z_pre_test_opti)*10^-9))
dlmwrite('pristine2_post_test_opti_crop.txt', ((Z_post_test_centered)*10^-9))

rs = 50


% center_4 = [186,465];
center_4 = [482,449];
range_2 = 30; % select how large you want the region to be: + and - the range in both axis

Z_diff_asperity = Z_diff(center_4(1)-range_2:center_4(1)+range_2,center_4(2)-range_2:center_4(2)+range_2);
Z_pre_test_opti_asp_iso = Z_pre_test_opti(center_4(1)-range_2:center_4(1)+range_2,center_4(2)-range_2:center_4(2)+range_2);
Z_post_test_centeredasp_iso = Z_post_test_centered(center_4(1)-range_2:center_4(1)+range_2,center_4(2)-range_2:center_4(2)+range_2);

center_5 = [31, 29]
range_3 = 6
Z_diff_asperity_loc = Z_diff_asperity(center_5(1)-range_3:center_5(1)+range_3,center_5(2)-range_3:center_5(2)+range_3);
Z_diff_asperity_loc_reshaped = reshape(Z_diff_asperity_loc, length(Z_diff_asperity_loc(:,1))*length(Z_diff_asperity_loc(1,:)),1);
Volume_Z_diff_asperity_loc = sum(Z_diff_asperity_loc_reshaped)

cross_section_before = Z_pre_test_opti_asp_iso(33,:)-65;
cross_section_after = Z_post_test_centeredasp_iso(35,:)-45;
cross_section_length = 1:43.945:2680.645;

figure
image((Z_pre_test_opti_asp_iso-mean(mean(Z_pre_test_opti))+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((Z_post_test_centeredasp_iso-mean(mean(Z_post_test_centered))+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((Z_diff_asperity-mean(mean(Z_diff_asperity))+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()
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 2700 -30 35])
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 \it{L} (nm) ');  %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');


% 
% rs = 90;
% 
% figure
% image((Z_diff_abs-mean(mean(Z_diff_abs))+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

%% We define the new difference height map here and plot it again


Z_diff2 = Z_post_test_centered - Z_pre_test_opti

opti_diff2 = mean(mean(Z_diff2))


rs = 50;

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

%% Here the functions that we use fr the script


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 = pos2-pos1; % determine the position of id2 with respect to id1
% Normally pos1 = pos2
pos_1 = pos1;
pos_2 = pos2;
pos = pos_1;
% a = 550; % random interger, not used later
% b = 400; % same
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 = zeros (20);
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(h, g) = difference_tilt;
        tilt(itb) = difference_tilt;
%         
        itb = itb +1
% %         ita = ita +1
    end
    
end

out1 = tilt;

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
