%% demo for Hyperspectral Restoration and Fusion with Multispectral Imagery...
%%                 by Recasting Low-Rank Tensor Approximation

%--------------Brief description-------------------------------------------
%
% 
% This demo implements HSI restoration and  spatial-resolution improvement by LRTA [1]
%
%
% More details in:
%
% [1] N. Liu, L. Li, W. Li, R. Tao J. Fowler and J. Chanussot, "Hyperspectral Restoration and Fusion...
%                   with Multispectral Imagery by Recasting Low-Rank Tensor Approximation," 
%                   in IEEE Transactions on Geoscience and Remote Sensing
% http://dx.doi.org/10.1109/TGRS.2020.3049014

% also in [2] N. Liu, W. Li, R. Tao and J. Fowler, "Wavelet-Domain Low-Rank/Group-Sparse ...
% Destriping for Hyperspectral Imagery," in IEEE Transactions on Geoscience and Remote Sensing
% http://dx.doi.org/10.1109/TGRS.2019.2933555
%
%
% If there are any problems to use the code and if you need more details about the data, code and paper, 
%  please feel free to contact us.
%
% contact: liwei089@ieee.org (Wei Li),liuna@ieee.org (Na Liu),1141017448@qq.com (Lu Li)

clear all;clc;close all;
addpath(genpath('D:\00STUDY\03_Fusion_LRTA\demo-LRTA-v2'))
addpath(genpath('Data'))
addpath(genpath('Methods'))
load('PaviaU.mat');
load('R.mat');


S = paviaU(1:256,1:256,11:end);
S = S./max(S(:));
[M N L] = size(S);

ratio = 8;
sz = [M N];
s0 = 1;
  psf        =    fspecial('gaussian',ratio,2);
  par.fft_B      =    psf2otf(psf,sz);
  par.fft_BT     =    conj(par.fft_B);
  par.H          =    @(z)H_z(z, par.fft_B, ratio, sz,s0 );
  par.HT         =    @(y)HT_y(y, par.fft_BT, ratio, sz,s0);

  F=R;
F=F(:,1:end-10);
for band = 1:size(F,1)
        div = sum(F(band,:));
        for i = 1:size(F,2)
            F(band,i) = F(band,i)/div;
        end
end
  
S_bar = hyperConvert2D(S);
hyper= par.H(S_bar);
MSI = hyperConvert3D((F*S_bar), M, N);

%
D=spatial_samples(psf,M,N)';
S3=Unfold(S,[M,N,L],3);
HSI_3=S3*D;
HSI=Fold(HSI_3,[M/ratio,N/ratio,L],3); 
  
%% degradation add stripes  
[H,W,Dim]=size(S);

%% add stripes
W1 = ones(H/size(psf,1),W/size(psf,2),Dim);
W2=ones(H,W,4);

LSHS = W1.*HSI;
MSI = W2.*MSI;

rate = 0.6;
Mean=0.2;%std(DataTest(:))% (10/255)/max(DataTest(:));%;
[Is, stripe_locations]  =  make_stripes(LSHS,rate,Mean);

 W1(:,stripe_locations) = 0;

%% LRTA with completion

N_max=max(M,max(N,L));
  alphia(1)=sqrt(M/N_max);
  alphia(2)=sqrt(N/N_max);
  alphia(3)=100*sqrt(L/N_max);
  alphia = alphia / sum(alphia);
  

 mu=0.01; 
 beta=0.5;
 gamma=0.5;
  maxIter=60;
 t0=clock;

Z_LRTA = RS_img_Fusion_mask_Fast(MSI,Is,D,F,W1,W2,alphia,mu,beta,gamma,maxIter,ratio);
t_LRTA=etime(clock,t0);
[psnr_LRTA,rmse_LRTA, ergas_LRTA, sam_LRTA, uiqi_LRTA,ssim_LRTA,DD_LRTA,CC_LRTA] = quality_assessment(double(im2uint8(S)), double(im2uint8(Z_LRTA)), 0, 1.0/ratio);
[psnr_LRTA_stripe,rmse_LRTA_stripe, ergas_LRTA_stripe, sam_LRTA_stripe, uiqi_LRTA_stripe,ssim_LRTA_stripe,DD_LRTA_stripe,CC_LRTA_stripe] = ...
    quality_assessment(double(im2uint8(HSI)), double(im2uint8(Is)), 0, 1.0/ratio);





