%% The code is written by Behnaz Pirzamanbein, bepi@dtu.dk last version 2020.08.04
%% measuring the distance betwee two sruface and length of each surface
addpath(genpath('.\Functions'))

save_folder = 'Results\measure\';
if ~exist(save_folder,'dir')
    mkdir(save_folder)
end 


%% load the surfaces
samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'};
degree = {'0','45','90','180'};

% if you want to run the measure separately you can use the for loop
% for i = 1:4
%     for j = 1:4
        sample = samples{i};
        dg = degree{j};

        fprintf('%s_%s \n',sample,dg)
        
        inner = load(['Results\layer\inner_converted_layer_',sample,'_',dg]);
        outer = load(['Results\layer\outer_converted_layer_',sample,'_',dg]);       

        %% smooth the layers
        disp('Smoothing Layers ...')
        
        p1 = mean(inner.surface_XYZ,4);
        p2 = mean(outer.surface_XYZ,4);

        Ba = regularization_matrix_open(size(p1,1),100,5000);
        p1_smooth = [];
        p2_smooth = [];
        for m = 1:size(p1,2)
            p1_smooth(:,m,:) = (squeeze(p1(:,m,:))'*Ba')';
            p2_smooth(:,m,:) = (squeeze(p2(:,m,:))'*Ba')';
        end

        Bb = regularization_matrix_open(size(p1,2),1000,5000);
        p1_s = [];
        p2_s = [];
        for n = 1:size(p1,1)
            p1_s(n,:,:) = (squeeze(p1_smooth(n,:,:))'*Bb')';
            p2_s(n,:,:) = (squeeze(p2_smooth(n,:,:))'*Bb')';
        end
        
        clear p1_smooth p2_smooth

        %% Euclidean distances
        disp('Computing Euclidean Distances ...')

        dist = [];
        for s = 1:size(p1,2)
            dist(:,s) = sqrt(sum((p1_s(:,s,1:2) - p2_s(:,s,1:2)).^2,3));
        end

        if flag_save
            save([save_folder,'Dist_',sample,'_',dg],'dist')
        end

%     end
% end
close all