function [initial_tfm_matrix]=transformed_matrix_initial(args)
% -------------------------output------------------------------------------
% tfm_matrix:   resulted transform matrix.

original_args=args;                                                               % back up

%% cut the image
expand_rate=0.8;                                                                   % empirical parameter 
initial_points=args.initial_points;
left_bound=ceil(max(initial_points(1, 1)-expand_rate*(initial_points(1, 2)-initial_points(1, 1)), 1));
right_bound=floor(min(initial_points(1, 2)+expand_rate*(initial_points(1, 2)-initial_points(1, 1)), size(args.input_image, 2)));
top_bound=ceil(max(initial_points(2, 1)-expand_rate*(initial_points(2, 2)-initial_points(2, 1)), 1));
bottom_bound=floor(min(initial_points(2, 2)+expand_rate*(initial_points(2, 2)-initial_points(2, 1)), size(args.input_image, 1)));
new_image=zeros(bottom_bound-top_bound+1, right_bound-left_bound+1,  size(args.input_image, 3));                
for c=1:size(args.input_image, 3)
    new_image(:, :, c)=args.input_image(top_bound:bottom_bound, left_bound:right_bound, c);                         % for aviod boundary effect
end
args.input_image=uint8(new_image);
args.center=args.center+[1-left_bound; 1-top_bound];                                                                % obtain the center of boundary extension image
%% down-sample the image if the focus is too large.
pre_scale_matrix=eye(3);
focus_threshold=200; % The max length of the focus is 100, if too large downsample to matrix of this size.
min_length=min(original_args.focus_size);
if min_length>focus_threshold
    s=min_length/focus_threshold;
    dst_pt_size=round([bottom_bound-top_bound+1, right_bound-left_bound+1]/s);
%     args.input_image=imfilter(args.input_image, fspecial('gaussian', ceil(2^(s-1)), ceil(2^(s-1))));
    args.input_image=imresize(args.input_image, dst_pt_size);
    pre_scale_matrix=diag([s s 1]);
end

%% adjust the parsed parameters.
initial_tfm_matrix=args.initial_tfm_matrix;
initial_tfm_matrix=inv(pre_scale_matrix)*initial_tfm_matrix*pre_scale_matrix;
args.initial_tfm_matrix=initial_tfm_matrix;
args.focus_size=round(args.focus_size/pre_scale_matrix(1, 1));
args.center=args.center/pre_scale_matrix(1, 1);

if args.branch==1
    %% Branch-and-bound if set
    %% step 1: prepare data for the lowest resolution.
    total_scale=floor(log2(min(args.focus_size)/args.focus_threshold));
    downsample_matrix=[0.5 0 0; 0 0.5 0; 0 0 1];
    scale_matrix=downsample_matrix^total_scale;
    tfm=maketform('projective', scale_matrix');
    if args.blur
       %% blur if set
        input_image=imfilter(args.input_image, fspecial('gaussian', ceil(args.blur_kernel_size_k*2^total_scale), ceil(args.blur_kernel_sigma_k*2^total_scale)));
    else
        input_image=args.input_image;
    end
    input_image=imtransform(input_image, tfm, 'bicubic');                                                     % bicubic zoom the input(why?)
    if size(input_image, 3)>1
        input_image=rgb2gray(input_image);
    end
    input_image=double(input_image);
    initial_tfm_matrix=scale_matrix*args.initial_tfm_matrix*inv(scale_matrix);
    center=floor(transform_point(args.center, scale_matrix));                                                 % the center of the bicubic zoomed boundary extension image(dismatching with the next focus center?)
    focus_size=floor(args.focus_size/2^total_scale);                                                          % focus size of the bicubic zoomed original image(original image*scale)
    f_branch=zeros(3, 2*args.branch_accuracy+1);
    Dotau_branch=cell(3, 2*args.branch_accuracy+1);
    result_tfm_matrix=cell(3, 2*args.branch_accuracy+1);
    %% step 2: design branch-and-bound method.
            max_rotation=args.branch_max_rotation;
            max_skew=args.branch_max_skew;
%             level=3;
%             candidate_matrix=cell(3, 2*args.branch_accuracy+1);
            level=1;
            candidate_matrix=cell(1, 2*args.branch_accuracy+1);
            for i=1:2*args.branch_accuracy+1
                candidate_matrix{1, i}=eye(3);
                theta=-max_rotation+(i-1)*max_rotation/args.branch_accuracy;
                candidate_matrix{1, i}(1:2, 1:2)=[cos(theta) -sin(theta); sin(theta) cos(theta)];
%                 candidate_matrix{2, i}=eye(3);
%                 candidate_matrix{2, i}(1, 2)=-max_skew+(i-1)*max_skew/args.branch_accuracy;
%                 candidate_matrix{3, i}=eye(3);
%                 candidate_matrix{3, i}(2, 1)=-max_skew+(i-1)*max_skew/args.branch_accuracy;
            end

    gap=5;
    BLACK_MATRIX=zeros(focus_size(1)*level+gap*(level-1), focus_size(2)*(2*args.branch_accuracy+1)+gap*2*args.branch_accuracy);
    %% step 3: begin branch-and-bound
    normal_outer_max_iter=args.outer_max_iter;
    normal_display_inter=args.display_result;
    args.outer_max_iter=1; % for debug, set it to 1;
    args.display_result=0;
    for i=1:level
        for j=1:2*args.branch_accuracy+1
           tfm_matrix=inv(candidate_matrix{i, j}*inv(initial_tfm_matrix));
           args.figure_no=(i-1)*level+j;
           args.save_path=[];

           %% update 1: assume that there's no great error, we directly
           %% compute nuclear norm of Dotau as the selection criterion record result
            image_size=size(input_image);
            image_center=floor(center);
            focus_center=zeros(2, 1);
            focus_center(1)=floor((1+focus_size(2))/2);
            focus_center(2)=floor((1+focus_size(1))/2);
            UData=[1-image_center(1) image_size(2)-image_center(1)];
            VData=[1-image_center(2) image_size(1)-image_center(2)];
            XData=[1-focus_center(1) focus_size(2)-focus_center(1)];
            YData=[1-focus_center(2) focus_size(1)-focus_center(2)];
            tfm=fliptform(maketform('projective', tfm_matrix'));
            Dotau=imtransform(input_image, tfm, 'bilinear', 'XData', XData, 'YData', YData, 'UData', UData, 'VData', VData, 'Size', focus_size);
            Dotau=Dotau/norm(Dotau, 'fro');
            [U S V]=svd(Dotau);
            f=sum(sum(S));
%            disp(['branching: level=', num2str(i), ', idx=', num2str(j)]);
           start=[(focus_size(1)+gap)*(i-1)+1, (focus_size(2)+gap)*(j-1)+1];
           BLACK_MATRIX(start(1):(start(1)+focus_size(1)-1), start(2):(start(2)+focus_size(2)-1))=Dotau;
           f_branch(i, j)=f;
           Dotau_branch{i, j}=Dotau;
           result_tfm_matrix{i, j}=tfm_matrix;
        end
        [value index]=min(f_branch(i, :));
        initial_tfm_matrix=result_tfm_matrix{i, index};
    end
    %% step 4: adapt initial_tfm_matrix to highest-resolution.
    initial_tfm_matrix=inv(scale_matrix)*initial_tfm_matrix*scale_matrix;
    
    %% step 5: show inter result if necessary.
    if normal_display_inter==1
        figure; imshow(BLACK_MATRIX, [], 'DisplayRange', [0 max(max(max(BLACK_MATRIX)))]);
    end
    position_matrix =[(focus_size(2)+gap)*(index-1)+1,(focus_size(2)+gap)*(index)-gap;1,focus_size(1)];  
    X1=[position_matrix(1,1) position_matrix(1,1) position_matrix(1,2) position_matrix(1,2) position_matrix(1,1)];
    Y1=[position_matrix(2,1) position_matrix(2,2) position_matrix(2,2) position_matrix(2,1) position_matrix(2,1)];
    hold on;
    plot(X1, Y1, 'r-', 'LineWidth',2);title('Rotation estimation for initialization')
    args.outer_max_iter=normal_outer_max_iter;
    args.display_result=normal_display_inter;
end