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

%--------------Brief description-------------------------------------------
%
% 
% This demo implements HSI 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
%
%
% 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);
HSI = hyperConvert3D(hyper,M/ratio, N/ratio );


%% LRTA
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);

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_Fast(MSI,HSI,D,alphia,mu,beta,gamma,F,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);

