function [T,S]=RS_img_Fusion_mask(T0,HSMS,LSHS,PSF,S,W1,W2,alphia,mu,beta,gamma)
%%W1--mask of HSMS, W2--mask of LSHS
addpath('D:\14_TC\TC_fusion\PROPACK');
[H0,W0,diml]=size(HSMS);
[Hl,Wl,Dim]=size(LSHS);
[TH,TW]=size(PSF);
norm_T0_fro=norm(T0(:),'fro');
%%%%%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
    padding_t=sum(sum(LSHS(:,:,i).*W1(:,:,i)))/nnz(W1(:,:,i));
    temp= squeeze(LSHS(:,:,i));
    temp(squeeze(W1(:,:,i))==0)=padding_t;
    LSHS(:,:,i)=temp;
end
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');

%%%%% initalzation S
flag=1;
if isempty(S)
    flag=0;
    S=zeros(Dim,diml);
    for i=1:diml
        S(:,i)=lsqnonneg(M3,HSMS_3(:,i));
    end%M3\HSMS_3;
end
%%%%% 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;%500;
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));

W1_3=Tensor_unfold(W1,3);
W2_3=Tensor_unfold(W2,3);

%% Start main loop
convergenced = 0;
iter = 0;
while iter<maxIter
    iter=iter+1;
    Tk=T;
    %%update M1
    
    %[M1,sv1]=svd_threshold(T1+Y1/mu,sv1,alphia(1)/mu,opt);
    [M1, ~, ~] = Pro2TraceNorm(T1+Y1/mu,alphia(1)/mu);
    %%update M2
    T2=Unfold(T,[H,W,Dim],2);
    [M2, ~, ~]=Pro2TraceNorm(T2+Y2/mu,alphia(2)/mu);
    %%update M3
    T3=Unfold(T,[H,W,Dim],3);
    delta_f_M3=-mu*(T3-M3+Y3/mu)+beta*(M3*DTD-DTL+Y_1*D'/beta)+gamma*(STS*M3-HST+SRF'*Y_2/gamma);
    [M3, ~, ~]=Pro2TraceNorm(M3-delta_f_M3/tau,alphia(3)/tau);
    
    
    
    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*D'*(W1_3.*(D*M3-LSHS_3+Y_1/beta))+gamma*(W2_3.*(M3*S-HSMS_3+Y_2/gamma))*S';
    [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
    if ~flag
        for i=1:diml
            S(:,i)=lsqnonneg(M3,HSMS_3(:,i)-Y_2(:,i)/gamma);
        end
    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*(W1_3.*(D*M3-LSHS_3));
    Y_2=Y_2+gamma*(W2_3.*(M3*S-HSMS_3));
    %%%% 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
        TruthErrrel=T-T0;
        disp(['iter ' num2str(iter) ',TruthErrrel=' num2str(norm(TruthErrrel(:),'fro')/norm_T0_fro) ',relChT=' num2str(relChgT)...
            ',recErr=' num2str(recErr)]);
    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