PrepRegistration.m

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