function [I_patch, R_Patch,  tfm_matrix]=RTLR(input_image,initial_tfm_matrix, para)

% -------------------------input------------------------------------------
% input_image:   extract image patch;
% initial_rain:  extract rain patch;
% initial_tfm_matrix:   3-by-3 matrix, optional, the initial transform
%                       matrix. If not specified set it to be eye(3);
% para:         parameters for both inner-loop and outer-loop, must be specified!
% -------------------------output-----------------------------------------
% I_patch:            estimated clear background
% R_Patch:            estimated rain streak
% Err_patch:          sparse error
% tfm_matrix:         resulted transform matrix.


%% parse parameters.

input_image=double(input_image);
center     = para.center;
focus_size = para.focus_size;

if isempty(para.outer_tol)
    outer_tol=5e-5;
else
    outer_tol=para.outer_tol;
end

if isempty(para.outer_max_iter)
    outer_max_iter=50;
else
    outer_max_iter=para.outer_max_iter;
end

inner_para=[];
if isempty(para.inner_tol)
    inner_para.tol=1e-7;
else
    inner_para.tol=para.inner_tol;
end

if isempty(para.inner_c)
    inner_para.c=1;
else
    inner_para.c=para.inner_c;
end

if isempty(para.inner_alpha)
    inner_para.alpha=[];
else
    inner_para.alpha=para.inner_alpha;
end

if isempty(para.inner_display_period)
    inner_para.display_period=100;
else
    inner_para.display_period=para.inner_display_period;
end

if isempty(para.inner_max_iter)
    inner_para.max_iter=inf;
else
    inner_para.max_iter=para.inner_max_iter;
end

%% decide origin of the two axes.
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)];

Dotau_series={};

%% prepare initial data
input_du=imfilter(input_image, -fspecial('sobel')'/8);
input_dv=imfilter(input_image, -fspecial('sobel')/8);
tfm_matrix=initial_tfm_matrix;
tfm=fliptform(maketform('projective', tfm_matrix'));
Dotau=imtransform(input_image, tfm, 'bilinear', 'UData', UData, 'VData', VData, 'XData', XData, 'YData', YData, 'size', focus_size);

Dotau_series{1}=Dotau;
initial_image=Dotau;
du=imtransform(input_du, tfm, 'bilinear', 'UData', UData, 'VData', VData, 'XData', XData, 'YData', YData, 'size', focus_size);
dv=imtransform(input_dv, tfm, 'bilinear', 'UData', UData, 'VData', VData, 'XData', XData, 'YData', YData, 'size', focus_size);
% du= du/norm(Dotau, 'fro')-(sum(sum(Dotau.*du)))/(norm(Dotau, 'fro')^3)*Dotau;
% dv= dv/norm(Dotau, 'fro')-(sum(sum(Dotau.*dv)))/(norm(Dotau, 'fro')^3)*Dotau;
% Dotau=Dotau/norm(Dotau, 'fro');
du= du/255;
dv= dv/255;
Dotau=Dotau/255;
tau=tfm2para(tfm_matrix, XData, YData, 'affine');
J=jacobi(du, dv, XData, YData, tau, 'affine');
S=constraints(tau, XData, YData, 'affine');
outer_round=0;
%% begin main loop.
while 1
   outer_round=outer_round+1;
   tic
   [I_patch, R_Patch, delta_tau] = ADMM_RTLR2(Dotau, J, S, inner_para);
%    [I_patch, R_Patch, delta_tau] = ADMM_RTLR3(Dotau, J, S, inner_para);
%    [I_patch, R_Patch, Err_Patch, delta_tau] = ADMM_RTLR1(Dotau, J, S, inner_para);
   toc
% figure,imshow(Dotau,[])
% figure,imshow(Err_Patch,[])
% figure,imshow(I_patch,[])
% figure,imshow(R_Patch,[])
   %% update Dotau.
   delta_tau =0;
   tau=tau+delta_tau;
   tfm_matrix=para2tfm(tau, XData, YData, 'affine');
   tfm=fliptform(maketform('projective', tfm_matrix'));
   Dotau=imtransform(input_image, tfm, 'bilinear', 'UData', UData, 'VData', VData, 'XData', XData, 'YData', YData, 'size', focus_size);
   Dotau_series{outer_round+1}=Dotau;
   %% judge convergence
   if outer_round>=outer_max_iter
       break;
   end
   %% record data and prepare for the next round.
   du=imtransform(input_du, tfm, 'bilinear', 'UData', UData, 'VData', VData, 'XData', XData, 'YData', YData, 'size', focus_size);
   dv=imtransform(input_dv, tfm, 'bilinear', 'UData', UData, 'VData', VData, 'XData', XData, 'YData', YData, 'size', focus_size);
%    du= du/norm(Dotau, 'fro')-(sum(sum(Dotau.*du)))/(norm(Dotau, 'fro')^3)*Dotau;
%    dv= dv/norm(Dotau, 'fro')-(sum(sum(Dotau.*dv)))/(norm(Dotau, 'fro')^3)*Dotau;
%    Dotau=Dotau/norm(Dotau, 'fro');
   du= du/255;
   dv= dv/255;
   Dotau=Dotau/255;
   J=jacobi(du, dv, XData, YData, tau, 'affine');
   S=constraints(tau, XData, YData, 'affine');
end
