%% The code is written by Behnaz Pirzamanbein, bepi@dtu.dk last version 2020.08.04
%% summing over 10 pixels along normal direction of the smoothed surface
addpath(genpath('.\Functions'))

save_folder = 'Results\registration\';
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};

            inner = load(['Results\layer\inner_converted_layer_',sample,'_',dg]);
            outer = load(['Results\layer\outer_converted_layer_',sample,'_',dg]);

            flag_layer = 'inner';
            load(['Results\point\',flag_layer,'_points_',sample,'_',dg])

            %% Contrast enhansed volume
            switch answer
                case 'No'
                    Vout = Vol;
                case 'Yes'
                    ce = points.contrast_enhansed_range;
                    Vout = (Vol-ce(1))/(ce(2)-ce(1))*ce(2); 
                    Vout(Vout < ce(1)) = ce(1); 
                    Vout(Vout > ce(2)) = ce(2);
            end
            %% smooth the 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

            %% sum up over 10 step in normal direction inward of the material
            no_slice = points.slice_no(end)-points.slice_no(1)+1;

            s1_conv = [];
            s2_conv = [];

            h = waitbar(0,'1','Name','Extracting layer ...');
            for s = 1:size(p1,2)
                slice_no = points.slice_no(1)+s-1;

                % inner surface
                range = 5:-1:-10;
                normals = compute_normals(p1_s(:,s,1:2),1);

                % grid along normals
                p_X = p1_s(:,s,1) + normals(:,1)*range;
                p_Y = p1_s(:,s,2) + normals(:,2)*range;

                interp_method = 'nearest';
                s1_conv(:,:,s) = interp2(Vout(:,:,slice_no),p_X,p_Y,interp_method);


                % outer surface
                range = 10:-1:-10;
                normals = compute_normals(p2_s(:,s,1:2),1);

                % grid along normals
                p_X = p2_s(:,s,1) + normals(:,1)*range;
                p_Y = p2_s(:,s,2) + normals(:,2)*range;


                interp_method = 'nearest';
                s2_conv(:,:,s) = interp2(Vout(:,:,slice_no),p_X,p_Y,interp_method);

                waitbar(s/no_slice,h,sprintf('slice: %d',slice_no))
            end
            close(h)

            s1 = squeeze(sum(s1_conv,2));
            s2 = squeeze(sum(s2_conv,2));

            if flag_save
                save([save_folder,'registration_',sample,'_',dg],'s1','s2')
                imwrite(rescale(s1),[save_folder,'registration_s1_',sample,'_',dg,'.png'])
                imwrite(rescale(s2),[save_folder,'registration_s2_',sample,'_',dg,'.png'])
            end
%     end
% end

close all