function [X]   =   ImUpdata(Im, X, FPN, Reconstruction, Weight, opts )
%% initialization
[m,n] = size(Im);
otfDx = psf2otf([1 -1],[m,n]);
otfDy = psf2otf([1;-1],[m,n]);

Dx  = zeros(m,n);
Dy  = zeros(m,n);
Z     = zeros(m,n);
Bx  = zeros(m,n);
By  = zeros(m,n);
Bz    = zeros(m,n);
[conjoDx,conjoDy,Denom1,Denom2] = getC(Im);

%% out-loop iteration
for iter = 1:opts.Innerloop_X
    %% X - subproblem
    aa = Reconstruction./(Weight+eps);    
    Denom  =   opts.belta_1 * Denom1 + opts.belta_2 * Denom2 + opts.gamma + opts.kappa;  
    Fx     =   ( opts.kappa * fft2( aa - FPN ) + opts.belta_1 * conjoDx.*fft2(Dx - Bx) + opts.belta_2 * conjoDy.*fft2(Dy - By) + opts.gamma * fft2(Z - Bz) + opts.belta_2  * Denom2.*fft2(aa))./Denom;     %% destriping 
    X      =   real( ifft2( Fx ) );  

    %% Z - subproblem
    Z      =   (opts.rho*Reconstruction + opts.gamma * X + opts.gamma * Bz)./(opts.gamma + opts.rho*Weight);              %% denoising

    %% Dx_s  and Dy_p - subproblem
     
    Du_y   =   real(ifft2(otfDy.*fft2(X - aa)));                           % row diff anisotropic TV
    Dy    =   soft_L1_shrink(By + Du_y, opts.lambda_2/opts.belta_2);
    
    Du_x   =   real(ifft2(otfDx.*fft2(X)));                                % column diff
    Dx     =   soft_L1_shrink(Bx + Du_x, opts.lambda_1/opts.belta_1);    

    %% ¸üÐÂ Bx, By, Bz
    
    Bx     =    Bx + Du_x - Dx;
    By     =    By + Du_y - Dy;  
    Bz     =    Bz + X - Z;
    
end

function [conjoDx,conjoDy,Denom1,Denom2] = getC(f)
% compute fixed quantities
sizeI = size(f);
otfDx = psf2otf([1,-1],sizeI);
otfDy = psf2otf([1;-1],sizeI);
conjoDx = conj(otfDx);
conjoDy = conj(otfDy);
Denom1 = abs(otfDx).^2;
Denom2 = abs(otfDy).^2;