%=====================================================================
% Programmer: Chih-Hsiang Hsieh
% E-mail: s104064515@m104.nthu.edu.tw
% Date: 2017/10/19
% -------------------------------------------------------
% Reference:
% C.-H. Lin, F. Ma, C.-Y. Chi, and C.-H. Hsieh,
% ``A convex optimization based coupled non-negative matrix factorization algorithm for hyperspectral and multispectral data fusion,"
% accepted by IEEE Trans. Geoscience and Remote Sensing, 2017.
%======================================================================
% A convex optimization based coupled NMF algorithm for hyperspectral superresolution via big data fusion
% [Z_fused,time] = ConvOptiCNMF(Yh,Ym,N,D,K)
%======================================================================
% Input
% Yh is low-spatial-resolution hyperspectral data cube of dimension rows_h*columns_h*M.
% (rows_h: vertical spatial dimension;
% columns_h: horizontal spatial dimension;
% M: spectral dimension.)
% Ym is high-spatial-resolution multispectral data cube of dimension rows_m*columns_m*Mm.
% (rows_m: vertical spatial dimension;
% columns_m: horizontal spatial dimension;
% M_m: spectral dimension.)
% N is the model order (should be greater than the number of endmembers).
% D is spectral response transform matrix of dimension M_m*M.
% K is blurring kernel of dimension k*k
% (The first dimension of K corresponds to the vertical spatial dimension;
% the second dimension of K corresponds to the horizontal spatial dimension.)
%----------------------------------------------------------------------
% Output
% Z_fused is super-resolved hyperspectral data cube of dimension rows_m*columns_m*M.
% time is the computational time (in secs).
%========================================================================
function [Z_fused,time] = ConvOptiCNMF(Yh,Ym,N,D,K)
t0=clock;
print_flag=0;
lambda1=0.001; lambda2=0.001; % SSD regularization parameter & L1-norm regularization parameter
eta=1; tildeeta=1; % penalty parameter in ADMM Algorithm 2 & penalty parameter in ADMM Algorithm 3
mode=1; % 0 uses Eq.(17)&(26) (for non-structured blur), 1 uses Lemmas 1 & 2 (for structured blur)
Max_iter=30; % maximum iteration of Algorithm 1
iterS=5; % maximum iteration of Algorithm 2
iterA=5; % maximum iteration of Algorithm 3
[rows_m,cols_m,bands_m]=size(Ym);
[rows_h,cols_h,M]=size(Yh);
Yh=reshape(Yh,[],M)';
[M,Lh]=size(Yh);
r1=rows_m/rows_h; % vertical dimension
r2=cols_m/cols_h; % horizontal dimension
r=max(ceil(r1),ceil(r2));
Ym=imresize(Ym,[r*rows_h,r*cols_h]);
Ym=reshape(Ym,[],bands_m)';
[Perm,g,B]=Permutation(K,r,rows_h,cols_h);
if mode==1, Ym=Ym*Perm'; end
I_M=speye(M);
P=sparse(zeros(0.5*M*N*(N-1),M*N));
x=1;
for n=1:N-1,
e_n=speye(N);
H_n=kron(e_n(:,n)',I_M);
for m=n+1:N,
P_m=kron(e_n(:,m)',I_M);
P(x:x+M-1,:)=H_n-P_m;
x=x+M;
end
end
PtP=P'*P; P=[];
[A,S_h,~]=HyperCSI_int(Yh,N); % initialization by HyperCSI [43]
S=S_h*kron(speye(Lh),ones(size(g')));
for k=1:Max_iter,
[S]=S_NonNeg_Lasso(N,Yh,Ym,A,S,D,g,B,iterS,lambda2,eta,mode); % Algorithm 2 (warm start)
[A]=A_ICE(N,Yh,Ym,S,A,D,g,B,iterA,lambda1,tildeeta,PtP,mode); % Algorithm 3 (warm start)
end
if mode==0,
Z_fused=A*S;
else
Z_fused=A*S*Perm;
end
Z_fused=reshape(Z_fused',r*rows_h,r*cols_h,M);
Z_fused=imresize(Z_fused,[rows_m,cols_m]); % resize back
time=etime(clock,t0);
%% subprogram 1 (implementation of Algorithm 2, solving the non-negative LASSO)
function [S]=S_NonNeg_Lasso(N,Yh,Ym,A,S_old,D,g,B,iterS,lambda2,eta,mode)
L=size(Ym,2);
y=reshape([reshape(Yh,[],1)',reshape(Ym,[],1)'],[],1);
x=reshape(S_old,[],1);
nu=sparse(zeros(N*L,1));
if mode==1,
C_bar1=kron(g',A)';
C_bar2=kron(speye(size(g,1)),D*A)';
C_bartC_bar=[C_bar1,C_bar2]*[C_bar1,C_bar2]';
C_bar_inv=C_bartC_bar+eta*speye(size(C_bartC_bar));
C_binv1=C_bar_inv\C_bar1;
C_binv2=C_bar_inv\C_bar2;
Ym_re=reshape(Ym,size(C_binv2,2),[]);
C_bYhYm=reshape(C_binv1*Yh+C_binv2*Ym_re,[],1);
for j=1:iterS,
% s update
xnu_re=reshape(eta*x-nu,size(C_bar_inv,2),[]);
s=reshape(C_bar_inv\xnu_re,[],1)+C_bYhYm;
% x update
x=s+nu/eta-lambda2/eta;
x(x<0)=0;
% nu update
nu=nu+eta*(s-x);
end
else % naive version
C=[kron(B',A)',kron(speye(L),D*A)']';
CtC=C'*C;
Cty=C'*y;
Lbs=chol(CtC+eta*speye(N*L),'lower');
for j=1:iterS,
% s update
q=Cty+eta*x-nu;
s=Lbs'\(Lbs\q);
% x update
x=s+nu/eta-lambda2/eta;
x(x<0)=0;
% nu update
nu=nu+eta*(s-x);
end
end
S=reshape(x,N,L);
return;
%% subprogram 2 (implementation of Algorithm 3, solving the ICE-regularized problem)
function [A]=A_ICE(N,Yh,Ym,S,A_old,D,g,B,iterA,lambda1,tildeeta,PtP,mode)
M=size(Yh,1);
y=reshape([reshape(Yh,[],1)',reshape(Ym,[],1)'],[],1);
z=reshape(A_old,[],1);
nu=zeros(M*N,1);
if mode==1,
S_3d=reshape(S,N,size(g,1),[]);
S_perm=permute(S_3d,[1,3,2]);
S_re=reshape(S_perm,[],size(g,1));
Sgv=reshape(S_re*g,N,[]);
CtCb1=kron(Sgv*Sgv',speye(M));
DtD=sparse(D'*D);
StS=S*S';
CtCb2=kron(StS,DtD);
CtC=CtCb1+CtCb2;
CtYh=reshape((Sgv*Yh')',[],1);
CtYm=reshape((S*(D'*Ym)')',[],1);
Cty=CtYh+CtYm;
CtCinv=CtC+lambda1*(PtP)+tildeeta*speye(M*N);
for j=1:iterA,
% a update
a=CtCinv\(Cty+tildeeta*z-nu);
% z update
z=a+nu/tildeeta;
z(z<0)=0;
% nu update
nu=nu+tildeeta*(a-z);
end
else % naive version
C=[kron((S*B)',sparse(eye(M)))' kron(S',D)']';
CtC=C'*C;
Cty=C'*y;
Lba=chol(CtC+lambda1*(PtP)+tildeeta*speye(MN),'lower');
for j=1:iterA,
% a update
q=Cty+tildeeta*z-nu;
a=Lba'\(Lba\q);
% z update
z=a+nu/tildeeta;
z(z<0)=0;
% nu update
nu=nu+tildeeta*(a-z);
end
end
A=reshape(z,M,N);
return;
%% subprogram 3
function [Pi,gv,B] = Permutation(K,r,rows_h,cols_h)
% This function helps user automatically generate the spatial spread
% transform matrix B from the blurring kernel K, automatically permute the
% pixels in a specific order (so that our fast closed-form solutions can be
% applied), and then accordingly revise the spatial spread transform matrix.
% This function considers the following 4 types of K, and automatically adapts
% it to the desired form with the original properties (variance) remained unchanged:
% 1. K is symmetric Gaussian, and k=sqrt((rows_m*columns_m)/(rows_h*columns_h))=r (the desired form).
% 2. K is symmetric Gaussian, but k does not equal to r.
% 3. K is not symmetric Gaussian but is uniform.
% 4. K is not symmetric Gaussian and is non-uniform.
blurkernel=size(K,1); r_square=r^2; [isgusyK,sigmaK]=isgusy(K);
if (all(K(:)==K(1))==1)
blurring_case=0;
elseif blurkernel==r && (isgusyK==1)
blurring_case=1;
elseif (isgusyK==1)
blurring_case=2;
else
blurring_case=3;
end
switch blurring_case
case 0
Temp=ones(r,r)/(r^2);
case 1
Temp=K;
case 2
Temp=fspecial('gaussian',[r,r],sigmaK);
case 3
zz=1;
ubound=150; lbound=0.05; interval=0.05;
a=zeros(((ubound-lbound)/interval)+1,1);
for sigma=lbound:interval:ubound
Temp=fspecial('gaussian',[blurkernel,blurkernel],sigma);
a(zz)=norm(Temp-K,'fro'); zz=zz+1;
end
[~,index]=min(a);
sigma=lbound+(index-1)*interval;
Temp=fspecial('gaussian',[r,r],sigma);
end
% generate B before permutation
Hv=speye(rows_h*cols_h);
for i=1:rows_h*cols_h
Hm=sparse(reshape(Hv(:,i),rows_h,cols_h));
B(:,i)=sparse(reshape(kron(Hm,Temp),[],1));
end
gv=reshape(Temp,[],1);
% start permut
