function [E_Img, Rain,  final_tfm_matrix]   =  RTLR_DeRaining( N_Img, ini_tfm_matrix, par )

    imacc = zeros(size(N_Img));
    imwgt = zeros(size(N_Img));
    [nrows,ncols]=size(N_Img);
    patchsize = par.patsize;
    pwidth = patchsize; 
    pheight = pwidth;
    rstep = floor(pheight/4);
    cstep = floor(pwidth/4);

    gfsize = double(int16(sqrt(2)*[pheight pwidth]/2)*2 + 5);
    gfilter1ch = fspecial('gaussian',[gfsize(1) gfsize(2)], mean(gfsize)/4);
    gfilter = gfilter1ch;
    steptimes = 2;

    for r = (0*rstep+int16(sqrt(2)*pheight/2)+5):steptimes*rstep:(nrows-int16(sqrt(2)*pheight/2)-5)
        for c = (0*cstep+int16(sqrt(2)*pwidth/2)+5):steptimes*cstep:(ncols-int16(sqrt(2)*pwidth/2)-5) 
            
          
            ridx = (r-floor(pheight/2)):(r+floor(pheight/2));
            cidx = (c-floor(pwidth/2)):(c+floor(pwidth/2));

 
            ridx2 = max(1, int16((r-sqrt(2)*pheight/2))):min(nrows, int16(r+sqrt(2)*pheight/2));
            cidx2 = max(1, int16((c-sqrt(2)*pwidth/2))):min(ncols, int16(c+sqrt(2)*pwidth/2));
            impatchlarge = N_Img(ridx2,cidx2);
            tfm=fliptform(maketform('projective', ini_tfm_matrix'));
            K=imtransform(impatchlarge, tfm, 'nearest');
            Dotau = K(ridx-r+ceil(size(K,1)/2),cidx-c+ceil(size(K,2)/2));
            
            input_du=imfilter(impatchlarge, -fspecial('sobel')'/8);
            input_dv=imfilter(impatchlarge, -fspecial('sobel')/8);
            du1=imtransform(input_du, tfm, 'nearest');
            dv1=imtransform(input_dv, tfm, 'nearest');
            du = du1(ridx-r+ceil(size(du1,1)/2),cidx-c+ceil(size(du1,2)/2));
            dv = dv1(ridx-r+ceil(size(dv1,1)/2),cidx-c+ceil(size(dv1,2)/2));
            focus_center(1)=floor((1+patchsize)/2);
            focus_center(2)=floor((1+patchsize)/2);
            XData=[1-focus_center(1) patchsize-focus_center(1)];
            YData=[1-focus_center(2) patchsize-focus_center(2)];
%             du= du/255;
%             dv= dv/255;
%             Dotau=Dotau/255;
            tau=tfm2para(ini_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]   =  LSSC_DeNoising( Dotau,J,S,par);
%                 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'));
                K=imtransform(impatchlarge, tfm, 'nearest');
                Dotau = K(ridx-r+ceil(size(K,1)/2),cidx-c+ceil(size(K,2)/2));
              %% judge convergence
                if outer_round>=par.outer_max_iter
                    break;
                end
                %% record data and prepare for the next round.
                du1=imtransform(input_du, tfm, 'nearest');
                dv1=imtransform(input_dv, tfm, 'nearest');
                du = du1(ridx-r+ceil(size(du1,1)/2),cidx-c+ceil(size(du1,2)/2));
                dv = dv1(ridx-r+ceil(size(dv1,1)/2),cidx-c+ceil(size(dv1,2)/2));
%                 du= du/255;
%                 dv= dv/255;
%                 Dotau=Dotau/255;
                J=jacobi(du, dv, XData, YData, tau, 'affine');
                S=constraints(tau, XData, YData, 'affine');
            end
           tfm=(maketform('projective', tfm_matrix'));
           impatchbw=imtransform(I_patch, tfm, 'nearest'); 
           [R_new,C_new ]= size(impatchbw);
           if mod(R_new,2) == 0 && mod(C_new,2) ~= 0
               impatchbw1 = zeros(R_new+1,C_new);
               impatchbw1(1:R_new,:)=impatchbw;
               impatchbw1(R_new+1,:)=impatchbw(R_new,:);
           elseif mod(C_new,2) == 0 && mod(R_new,2) ~= 0
               impatchbw1 = zeros(R_new,C_new+1);
               impatchbw1(:,1:C_new)=impatchbw;
               impatchbw1(:,C_new+1)=impatchbw(:,C_new);
           elseif mod(C_new,2) ==0 && mod(R_new,2) == 0
               impatchbw1 = zeros(R_new+1,C_new+1);
               impatchbw1(1:R_new,1:C_new)=impatchbw;
               impatchbw1(1:R_new,C_new+1)=impatchbw(:,C_new);
               impatchbw1(R_new+1,:) = impatchbw1(R_new,:);
           else
                impatchbw1 =impatchbw;
           end
           impatchbw = impatchbw1;
            hf = floor(size(impatchbw,1)/2);
            wf = floor(size(impatchbw,2)/2);
            gfilter2 = gfilter(int16((gfsize(1)+1)/2 - hf):int16((gfsize(1)+1)/2 + hf), int16((gfsize(2)+1)/2 - wf):int16((gfsize(2)+1)/2 + wf)); 
            gfilter2(impatchbw==0) = 0;
              
            imacc((r-hf):(r+hf),(c-wf):(c+wf)) = imacc((r-hf):(r+hf),(c-wf):(c+wf)) + gfilter2.*double(impatchbw);
            imwgt((r-hf):(r+hf),(c-wf):(c+wf)) = imwgt((r-hf):(r+hf),(c-wf):(c+wf)) + gfilter2;
        end
    end
    E_Img = round((imacc+1e-10)./(imwgt+1e-10)); 
%     E_Img = uint8(round((imacc+1e-10)./(imwgt+1e-10)));  