function T = RS_img_Fusion_mask_Fast(HSMS,LSHS,D,SRF,W1,W2,alphia,mu,beta,gamma,maxIter,ratio)
%   T=RS_img_Fusion_Fast(HSMS,LSHS,D,alphia,mu,beta,gamma,SRF,maxIter,ratio)
%%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);

%%%%% initalzation T by LSHS
T=imresize(LSHS,ratio);


norm_init_T=norm(T(:),'fro');
M1=Unfold(T,[H,W,Dim],1);
M2=Unfold(T,[H,W,Dim],2);
M3=Unfold(T,[H,W,Dim],3);

HSMS_3=Unfold(HSMS,[H,W,diml],3);%[H0,W0,diml]=size(HSMS);%Tensor_unfold(HSMS,3);
LSHS_3=Unfold(LSHS,[Hl,Wl,Dim],3);
norm_H=norm(HSMS_3,'fro');
norm_L=norm(LSHS_3,'fro');

%%%%% initalzation laplacian variable
Y1=zeros(H,W*Dim);
Y2=zeros(W,H*Dim);
Y3=zeros(Dim,H*W);
Y_1=zeros(Dim,Hl*Wl);
Y_2=zeros(diml,H*W);



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

s1=svds(D,1);
s2=svds(SRF,1);
tau=s1^2*beta+gamma*s2^2+mu;


% W1_3=Tensor_unfold(W1,3);
% W2_3=Tensor_unfold(W2,3);
W1_3=Unfold(W1,[Hl,Wl,Dim],3);
W2_3=Unfold(W2,[H0,W0,diml],3);
%% Start main loop
convergenced = 0;
iter = 0;
while iter<maxIter
    iter=iter+1;
    Tk=T;
    %%update M1
    T1=Unfold(T,[H,W,Dim],1);
    [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);
    delta_f_M3=-mu*(T3-M3+Y3/mu)+beta*(W1_3.*(M3*D-LSHS_3+Y_1/beta))*D'+gamma*SRF'*(W2_3.*(SRF*M3-HSMS_3+Y_2/gamma));
    % [M3,sv3]=svd_threshold(M3-delta_f_M3/tau,sv3,alphia(3)/tau,opt);
         
    [M3, ~, ~]=Pro2TraceNorm(M3-delta_f_M3/tau,alphia(3)/tau);
    %%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;   
    T=(Fold(M1-Y1/mu,[H,W,Dim],1)+Fold(M2-Y2/mu,[H,W,Dim],2)+Fold(M3-Y3/mu,[H,W,Dim],3))/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);
    
    Y1=Y1+mu*(Unfold(T,[H,W,Dim],1)-M1);
    Y2=Y2+mu*(Unfold(T,[H,W,Dim],2)-M2);
    Y3=Y3+mu*(Unfold(T,[H,W,Dim],3)-M3);   
    Y_1=Y_1+beta*(W1_3.*(M3*D-LSHS_3));
    Y_2=Y_2+gamma*(W2_3.*(SRF*M3-HSMS_3));
    
    


    
    %%%% convergence condition judgement
    diffT=T-Tk;
    relChgT=sqrt(diffT(:)'*diffT(:))/norm_init_T;
    T3=Unfold(T,[H,W,Dim],3);
    recErr=max(norm(T3*D-LSHS_3,'fro')/norm_L,norm(SRF*T3-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