function [T,RMS]=RS_img_Fusion(HSMS,LSHS,PSF,alphia,mu,beta,gamma,T0,S,maxIter)

[H0,W0,diml]=size(HSMS);
[Hl,Wl,Dim]=size(LSHS);
[TH,TW]=size(PSF);
%%%%%padding or cut HSMS
if H0<Hl*TH
    paddingH=Hl*TH-H0;
    HSMS=padarray(HSMS,[paddingH,0],'replicate','post');
else if H0>Hl*TH
        HSMS=HSMS(1:Hl*TH,:);
    end
end

if W0<Wl*TW
    paddingW=Wl*TW-W0;
    HSMS=padarray(HSMS,[0,paddingW],'replicate','post');
else if W0>Wl*TW
        HSMS=HSMS(:,1:Wl*TW);
    end
end
[H,W,~]=size(HSMS);
%%%%%initialzation
T=zeros(H,W,Dim);
PSF=PSF/sum(PSF(:));

D=spatial_samples(PSF,H,W);%%%%%computing sample matrix D by PSF

%%%%% initalzation T by LSHS
for i=1:Dim
    T(:,:,i)=kron(LSHS(:,:,i),PSF);
end
norm_init_T=norm(T(:),'fro');
M1=Tensor_unfold(T,1);
M2=Tensor_unfold(T,2);
M3=Tensor_unfold(T,3);
HSMS_3=Tensor_unfold(HSMS,3);
LSHS_3=Tensor_unfold(LSHS,3);
norm_H=norm(HSMS_3,'fro');
norm_L=norm(LSHS_3,'fro');
%S=zeros(Dim,diml);
%%%%% initalzation S
% for i=1:diml
%     S(:,i)=lsqnonneg(M3,HSMS_3(:,i));
% end%M3\HSMS_3;    
%%%%% initalzation laplacian variable
Y1=zeros(W*Dim,H);
Y2=zeros(H*Dim,W);
Y3=zeros(H*W,Dim);
Y_1=zeros(Hl*Wl,Dim);
Y_2=zeros(H*W,diml);

DTD=D'*D;
DTL=D'*LSHS_3;
STS=S*S';
HST=HSMS_3*S';

normfX = norm(M3,'fro');
tol1 = 1e-4;%threshold for the error in constraint
tol2 = 1e-5;%threshold for the change in the solutions
opt.tol = tol2;%precision for computing the partial SVD
%maxIter=100;
s1=svds(D,1);
s2=svds(S,1);
tau=s1^2*beta+gamma*s2^2+mu;

sv1=round(0.9*min(W*Dim,H));

sv2=round(0.9*min(H*Dim,W));

sv3=round(0.9*min(Dim,H*W));



%% Start main loop

% mu_max=5e-2;
% rho=1.1;
convergenced = 0;
iter = 0;
while iter<maxIter
    iter=iter+1;
    Tk=T;
    %%update M1
    T1=Tensor_unfold(T,1);
    [M1,sv1]=svd_threshold(T1+Y1/mu,sv1,alphia(1)/mu,opt);
    %%update M2
    T2=Tensor_unfold(T,2);
    [M2,sv2]=svd_threshold(T2+Y2/mu,sv2,alphia(2)/mu,opt);
    %%update M3
    T3=Tensor_unfold(T,3);
    delta_f_M3=-mu*(T3-M3+Y3/mu)+beta*(DTD*M3-DTL+D'*Y_1/beta)+gamma*(M3*STS-HST+Y_2*S'/gamma);
    [M3,sv3]=svd_threshold(M3-delta_f_M3/tau,sv3,alphia(3)/tau,opt);
    %%update T
    T=(fold_k(M1-Y1/mu,1,[H,W,Dim])+fold_k(M2-Y2/mu,2,[H,W,Dim])+fold_k(M3-Y3/mu,3,[H,W,Dim]))/3;    
    %%update S
%     for i=1:diml
%         S(:,i)=lsqnonneg(M3,HSMS_3(:,i)-Y_2(:,i)/gamma);
%     end
    %S=M3\(HSMS_3-Y_2/gamma);
    %%update laplacian variable
    Y1=Y1+mu*(Tensor_unfold(T,1)-M1);
    Y2=Y2+mu*(Tensor_unfold(T,2)-M2);
    Y3=Y3+mu*(Tensor_unfold(T,3)-M3);
    Y_1=Y_1+beta*(D*M3-LSHS_3);
    Y_2=Y_2+gamma*(M3*S-HSMS_3);
    
    %mu update
  %  mu=min(mu_max,rho*mu);
    
    %%%% convergence condition judgement
    diffT=T-Tk;
    relChgT=sqrt(diffT(:)'*diffT(:))/norm_init_T;
    M3=Tensor_unfold(T,3);
    recErr=max(norm(D*M3-LSHS_3,'fro')/norm_L,norm(M3*S-HSMS_3,'fro')/norm_H);
    convergenced = recErr <tol1 && relChgT < tol2;
    
    if iter==1 || mod(iter,50)==0 || convergenced
        disp(['iter ' num2str(iter) ',relChT=' num2str(relChgT)...
            ',recErr=' num2str(recErr)]);
        RMS=norm(T(:)-T0(:),'fro')/norm(T0(:),'fro')
      % RMS=0;
    end
    
    if convergenced
        break;
    end
end



end



function [out,sv]=svd_threshold(A,sv,thres,opt)
[U,S,V]=lansvd(A,sv, 'L', opt);
[H,W]=size(A);
n=min(H,W);
 S = diag(S);
svp = length(find(S>thres));
if svp < sv
    sv = min(svp + 1, n);
else
    sv = min(svp + round(0.05*n), n);
end

if svp>=1
    S = S(1:svp)-thres;
else
    svp = 1;
    S = 0;
end

U = U(:, 1:svp);
s = S;
V = V(:, 1:svp);

out = U*diag(s)*V';
end


function out=svd_threshold2(A,thres)
[U,S,V]=svd(A);
diag_S=max(diag(S)-thres,0);
out=U*diag(diag_S)*V';
end