%% The code is written by Behnaz Pirzamanbein, bepi@dtu.dk last version 2020.08.04
%% registration
samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'};

load_folder = 'Results\registration\';
save_folder = 'Results\registration\tform\';
if ~exist(save_folder,'dir')
    mkdir(save_folder)
end

% if you want to run the measure separately you can use the for loop
%for i = 1:4
    sample = samples{i};
    % loading data
    disp(sample)
    
    dg0_inner = imread([load_folder,'registration_s1_',sample,'_0.png']);
    dg45_inner = imread([load_folder,'registration_s1_',sample,'_45.png']);
    dg90_inner = imread([load_folder,'registration_s1_',sample,'_90.png']);
    dg180_inner = imread([load_folder,'registration_s1_',sample,'_180.png']);

    % estimate the resgitration T (using the generated function)
    Reg_dg0_inner = registerImages(dg0_inner,dg180_inner,30);
    Reg_dg45_inner = registerImages(dg45_inner,dg180_inner,30);
    Reg_dg90_inner = registerImages(dg90_inner,dg180_inner,30);
   
    save([save_folder,'reg_tfrom_inner_',sample],'*_inner')

    % loading data
    dg0_outer = imread([load_folder,'registration_s2_',sample,'_0.png']);
    dg45_outer = imread([load_folder,'registration_s2_',sample,'_45.png']);
    dg90_outer = imread([load_folder,'registration_s2_',sample,'_90.png']);
    dg180_outer = imread([load_folder,'registration_s2_',sample,'_180.png']);

    % estimate the resgitration T (using the generated function)
    Reg_dg0_outer = registerImages(dg0_outer,dg180_outer,30);
    Reg_dg45_outer = registerImages(dg45_outer,dg180_outer,30);
    Reg_dg90_outer = registerImages(dg90_outer,dg180_outer,30);  

    if flag_save
        save([save_folder,'reg_tfrom_outer_',sample],'*_outer')
    end
    
%end

%% check all against 180

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

load_folder = 'Results\registration\tform\';
save_folder = load_folder;
if ~exist(save_folder,'dir')
    mkdir(save_folder)
end

% if you want to run the measure separately you can use the for loop
%for i = 1:4
    sample = samples{i};
    inner = load([load_folder,'reg_tfrom_inner_',sample]);
    outer = load([load_folder,'reg_tfrom_outer_',sample]);

    MOVING = inner.dg0_inner;%imrotate(layers(:,:,1),180);
    movingRefObj = imref2d(size(inner.dg0_inner));

    fixedRefObj = inner.Reg_dg0_inner.SpatialRefObj;
    tform = inner.Reg_dg0_inner.Transformation;
    inner_0dg = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

    MOVING = inner.dg45_inner;%imrotate(layers(:,:,1),180);
    movingRefObj = imref2d(size(inner.dg45_inner));

    fixedRefObj = inner.Reg_dg45_inner.SpatialRefObj;
    tform = inner.Reg_dg45_inner.Transformation;
    inner_45dg = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

    MOVING = inner.dg90_inner;%imrotate(layers(:,:,1),180);
    movingRefObj = imref2d(size(inner.dg90_inner));

    fixedRefObj = inner.Reg_dg90_inner.SpatialRefObj;
    tform = inner.Reg_dg90_inner.Transformation;
    inner_90dg = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

    figure,
    subplot(221)
    imshowpair(inner_0dg,inner.dg180_inner),title('0 to 180')
    subplot(222)
    imshowpair(inner_45dg,inner.dg180_inner),title('45 to 180')
    subplot(223)
    imshowpair(inner_90dg,inner.dg180_inner),title('90 to 180')
    suptitle_modified(['inner using inner T ',sample],'none')
    if flag_save
        saveas(gcf,[save_folder,'inner_using_inner_T_',sample,'.fig'])
        saveas(gcf,[save_folder,'inner_using_inner_T_',sample,'.png'])
    end

    % check all outer agains 180 using inner transform T
    MOVING = outer.dg0_outer;%imrotate(layers(:,:,1),180);
    movingRefObj = imref2d(size(outer.dg0_outer));

    fixedRefObj = inner.Reg_dg0_inner.SpatialRefObj;
    tform = inner.Reg_dg0_inner.Transformation;
    outer_0dg = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

    MOVING = outer.dg45_outer;%imrotate(layers(:,:,1),180);
    movingRefObj = imref2d(size(outer.dg45_outer));

    fixedRefObj = inner.Reg_dg45_inner.SpatialRefObj;
    tform = inner.Reg_dg45_inner.Transformation;
    outer_45dg = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

    MOVING = outer.dg90_outer;%imrotate(layers(:,:,1),180);
    movingRefObj = imref2d(size(outer.dg90_outer));

    fixedRefObj = inner.Reg_dg90_inner.SpatialRefObj;
    tform = inner.Reg_dg90_inner.Transformation;
    outer_90dg = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

    figure,
    subplot(221)
    imshowpair(outer_0dg,outer.dg180_outer),title('0 to 180')
    subplot(222)
    imshowpair(outer_45dg,outer.dg180_outer),title('45 to 180')
    subplot(223)
    imshowpair(outer_90dg,outer.dg180_outer),title('90 to 180')
    suptitle_modified(['outer using inner T ',sample],'none')
    if flag_save
        saveas(gcf,[save_folder,'outer_using_inner_T_',sample,'.fig'])
        saveas(gcf,[save_folder,'outer_using_inner_T_',sample,'.png'])
    end
%end

close all

%% apply the inner T to the layers for measurement

samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'}; 

load_folder_Reg = 'Results\registration\tform\';
load_folder_layers = 'Results\layer\';
save_folder =  'Results\registration\layer\';
if ~exist(save_folder,'dir')
    mkdir(save_folder)
end

% if you want to run the measure separately you can use the for loop
%for i = 1:4
    clearvars -except samples i flag_data load_folder_Reg load_folder_layers save_folder
    close all
    
    sample = samples{i};
    fprintf('Load Data %s ...\n',sample)
    
    % for registrtaion
    inner = load([load_folder_Reg,'reg_tfrom_inner_',sample]);
    % for visualization
    outer = load([load_folder_Reg,'reg_tfrom_outer_',sample]);

    inner_layers_0 = load([load_folder_layers,'inner_layers_',sample,'_0']);
    outer_layers_0 = load([load_folder_layers,'outer_layers_',sample,'_0']);

    inner_layers_45 = load([load_folder_layers,'inner_layers_',sample,'_45']);
    outer_layers_45 = load([load_folder_layers,'outer_layers_',sample,'_45']);

    inner_layers_90 = load([load_folder_layers,'inner_layers_',sample,'_90']);
    outer_layers_90 = load([load_folder_layers,'outer_layers_',sample,'_90']);

    inner_layers_180 = load([load_folder_layers,'inner_layers_',sample,'_180']);
    outer_layers_180 = load([load_folder_layers,'outer_layers_',sample,'_180']);

    for s = 1:2
        fprintf('Registrtaion of Surface %d\n',s)
        
        dg0_inner = inner_layers_0.layers(:,:,s);
        dg45_inner = inner_layers_45.layers(:,:,s);
        dg90_inner = inner_layers_90.layers(:,:,s);
        dg180_inner = inner_layers_180.layers(:,:,s);

        dg0_outer = outer_layers_0.layers(:,:,s);
        dg45_outer = outer_layers_45.layers(:,:,s);
        dg90_outer = outer_layers_90.layers(:,:,s);
        dg180_outer = outer_layers_180.layers(:,:,s);

       
        MOVING = dg0_inner;
        movingRefObj = imref2d(size(dg0_inner));

        fixedRefObj = inner.Reg_dg0_inner.SpatialRefObj;
        tform = inner.Reg_dg0_inner.Transformation;
        [inner_0dg(:,:,s), inner_0dg_Ref] = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

        MOVING = dg45_inner;
        movingRefObj = imref2d(size(dg45_inner));

        fixedRefObj = inner.Reg_dg45_inner.SpatialRefObj;
        tform = inner.Reg_dg45_inner.Transformation;
        [inner_45dg(:,:,s), inner_45dg_Ref] = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

        MOVING = dg90_inner;
        movingRefObj = imref2d(size(dg90_inner));

        fixedRefObj = inner.Reg_dg90_inner.SpatialRefObj;
        tform = inner.Reg_dg90_inner.Transformation;
        [inner_90dg(:,:,s), inner_90dg_Ref] = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

        inner_180dg(:,:,s) = dg180_inner;

        MOVING = dg0_outer;
        movingRefObj = imref2d(size(dg0_outer));

        fixedRefObj = inner.Reg_dg0_inner.SpatialRefObj;
        tform = inner.Reg_dg0_inner.Transformation;
        [outer_0dg(:,:,s), outer_0dg_Ref] = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

        MOVING = dg45_outer;
        movingRefObj = imref2d(size(dg45_outer));

        fixedRefObj = inner.Reg_dg45_inner.SpatialRefObj;
        tform = inner.Reg_dg45_inner.Transformation;
        [outer_45dg(:,:,s), outer_45dg_Ref] = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

        MOVING = dg90_outer;
        movingRefObj = imref2d(size(dg90_outer));

        fixedRefObj = inner.Reg_dg90_inner.SpatialRefObj;
        tform = inner.Reg_dg90_inner.Transformation;
        [outer_90dg(:,:,s), outer_90dg_Ref] = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj);

        outer_180dg(:,:,s) = dg180_outer;

        figure,
        subplot(231)
        imshowpair(inner_0dg(:,:,s),inner.dg180_inner),title('0 to 180')
        ylabel('inner using inner T ','Interpreter','none')
        subplot(232)
        imshowpair(inner_45dg(:,:,s),inner.dg180_inner),title('45 to 180')
        subplot(233)
        imshowpair(inner_90dg(:,:,s),inner.dg180_inner),title('90 to 180')

        subplot(234)
        imshowpair(outer_0dg(:,:,s),outer.dg180_outer),title('0 to 180')
        ylabel('outer using inner T ','Interpreter','none')
        subplot(235)
        imshowpair(outer_45dg(:,:,s),outer.dg180_outer),title('45 to 180')
        subplot(236)
        imshowpair(outer_90dg(:,:,s),outer.dg180_outer),title('90 to 180')
        suptitle_modified(sample,'none')
        if flag_save
            saveas(gcf,sprintf('%sReg_s%d_%s.png',save_folder,s,sample))
        end
    end
    if flag_save
        disp('Save results ...')
        layers = inner_0dg;
        Ref = inner_0dg_Ref;
        save([save_folder,'Reg_inner_layers_',sample,'_0'],'layers','Ref')
        layers = inner_45dg;
        Ref = inner_45dg_Ref;
        save([save_folder,'Reg_inner_layers_',sample,'_45'],'layers','Ref')
        layers = inner_90dg;
        Ref = inner_90dg_Ref;
        save([save_folder,'Reg_inner_layers_',sample,'_90'],'layers','Ref')
        layers = inner_180dg;
        save([save_folder,'Reg_inner_layers_',sample,'_180'],'layers')

        layers = outer_0dg;
        Ref = outer_0dg_Ref;
        save([save_folder,'Reg_outer_layers_',sample,'_0'],'layers','Ref')
        layers = outer_45dg;
        Ref = outer_45dg_Ref;
        save([save_folder,'Reg_outer_layers_',sample,'_45'],'layers','Ref')
        layers = outer_90dg;
        Ref = outer_90dg_Ref;
        save([save_folder,'Reg_outer_layers_',sample,'_90'],'layers','Ref')
        layers = outer_180dg;
        save([save_folder,'Reg_outer_layers_',sample,'_180'],'layers')
    end
%end

close all 

%% apply the inner T to the mesh points and return the mesh back

samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'}; 

load_folder_Reg = 'Results\registration\tform\';
load_folder_layers = 'Results\layer\';
save_folder = 'Results\registration\layer\';
if ~exist(save_folder,'dir')
    mkdir(save_folder)
end

% if you want to run the measure separately you can use the for loop
%for i = 1:4
    clearvars -except samples i flag_data load_folder_Reg load_folder_layers save_folder
    close all
    
    sample = samples{i};
    fprintf('Load Data %s ...\n',sample)
    
    % for registrtaion
    inner = load([load_folder_Reg,'reg_tfrom_inner_',sample]);
    % for visualization
    outer = load([load_folder_Reg,'reg_tfrom_outer_',sample]);

    inner_layers_0 = load([load_folder_layers,'inner_layers_',sample,'_0']);
    outer_layers_0 = load([load_folder_layers,'outer_layers_',sample,'_0']);

    inner_layers_45 = load([load_folder_layers,'inner_layers_',sample,'_45']);
    outer_layers_45 = load([load_folder_layers,'outer_layers_',sample,'_45']);

    inner_layers_90 = load([load_folder_layers,'inner_layers_',sample,'_90']);
    outer_layers_90 = load([load_folder_layers,'outer_layers_',sample,'_90']);

    inner_layers_180 = load([load_folder_layers,'inner_layers_',sample,'_180']);
    outer_layers_180 = load([load_folder_layers,'outer_layers_',sample,'_180']);

    for s = 1:2
        % initial transformation (rotation or flip)
        fprintf('Registrtaion of the Mesh Points Surface %d\n',s)
        
        dg0_inner = inner_layers_0.layers(:,:,s);
        dg45_inner = inner_layers_45.layers(:,:,s);
        dg90_inner = inner_layers_90.layers(:,:,s);
        dg180_inner = inner_layers_180.layers(:,:,s);

        dg0_outer = outer_layers_0.layers(:,:,s);
        dg45_outer = outer_layers_45.layers(:,:,s);
        dg90_outer = outer_layers_90.layers(:,:,s);
        dg180_outer = outer_layers_180.layers(:,:,s);
        
        % 0 degree
        [tmp(:,2), tmp(:,1)] = find(dg0_inner);
        tform = inner.Reg_dg0_inner.Transformation;
        [inner_0dg(:,1,s), inner_0dg(:,2,s)] = transformPointsForward(tform, tmp(:,1), tmp(:,2));
        clear tmp Reg
                
        % 45 degree
        [tmp(:,2), tmp(:,1)] = find(dg45_inner);
        tform = inner.Reg_dg45_inner.Transformation;
        [inner_45dg(:,1,s), inner_45dg(:,2,s)] = transformPointsForward(tform, tmp(:,1), tmp(:,2));
        clear tmp Reg
        
        % 90 degree
        [tmp(:,2), tmp(:,1)] = find(dg90_inner);
        tform = inner.Reg_dg90_inner.Transformation;
        [inner_90dg(:,1,s), inner_90dg(:,2,s)] = transformPointsForward(tform, tmp(:,1), tmp(:,2));
        clear tmp Reg
        
        [inner_180dg(:,2,s), inner_180dg(:,1,s)] = find(dg180_inner);
        clear tmp Reg

        % 0 degree
        [tmp(:,2), tmp(:,1)] = find(dg0_outer);
        tform = inner.Reg_dg0_inner.Transformation;
        [outer_0dg(:,1,s), outer_0dg(:,2,s)] = transformPointsForward(tform, tmp(:,1), tmp(:,2));
        clear tmp Reg
        
        % 45 degree
        [tmp(:,2), tmp(:,1)] = find(dg45_outer);
        tform = inner.Reg_dg45_inner.Transformation;
        [outer_45dg(:,1,s), outer_45dg(:,2,s)] = transformPointsForward(tform, tmp(:,1), tmp(:,2));
        clear tmp Reg

        % 90 degree
        [tmp(:,2), tmp(:,1)] = find(dg90_outer);
        tform = inner.Reg_dg90_inner.Transformation;
        [outer_90dg(:,1,s), outer_90dg(:,2,s)] = transformPointsForward(tform, tmp(:,1), tmp(:,2));
        clear tmp Reg
        
        [outer_180dg(:,2,s), outer_180dg(:,1,s)] = find(dg180_outer);
        clear tmp Reg

    end
    if flag_save
        disp('Save results ...')
        points = inner_0dg;
        save([save_folder,'Reg_inner_points_',sample,'_0'],'points')
        points = inner_45dg;
        save([save_folder,'Reg_inner_points_',sample,'_45'],'points')
        points = inner_90dg;
        save([save_folder,'Reg_inner_points_',sample,'_90'],'points')
        points = inner_180dg;
        save([save_folder,'Reg_inner_points_',sample,'_180'],'points')

        points = outer_0dg;
        save([save_folder,'Reg_outer_points_',sample,'_0'],'points')
        points = outer_45dg;
        save([save_folder,'Reg_outer_points_',sample,'_45'],'points')
        points = outer_90dg;
        save([save_folder,'Reg_outer_points_',sample,'_90'],'points')
        points = outer_180dg;
        save([save_folder,'Reg_outer_points_',sample,'_180'],'points')
    end
%end

close all