-
Notifications
You must be signed in to change notification settings - Fork 2
/
plug_abundance.m
102 lines (93 loc) · 3.71 KB
/
plug_abundance.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
% Script for unmixing with Plug-and-Play Priors by denoising abundance maps
%
% Reference:
% Hyperspectral Unmixing Via Plug-And-Play Priors
% Xiuheng Wang, Min Zhao, Jie Chen
%
% A Plug-and-Play Priors Framework for Hyperspectral Unmixing
% Min Zhao, Xiuheng Wang, Jie Chen
%
% 2020/09/08
% Implemented by
% Min Zhao, Xiuheng Wang
clear;clc;
close all;
path(path,genpath(pwd));
recMode = { 'BM3D' 'WNNM' 'GSRC', 'AST-NLS', 'PGPD', ...
'MSEPLL', 'DnCNN', 'SSC_GSM' 'ACPT' , 'TWSC' ...
'NCSR' 'ACVA' 'GMM_EPLL' 'LMM_EPLL' 'GGMM_EPLL'};
recMode_id = 7;
ImgOrg = 1;
inPar.isShow = 0;
%% load data
Dict = importdata('real_data/end4_jas.mat');
a = Dict.A;
M = Dict.M;
y = importdata('real_data/jaspeRidge_nor.mat');
y = double(y);
%% Parameter settings
[L, N] = size(y);
p = size(M, 2);
SMALL = 1e-12;
I_1 = ones(p,1);
I_2 = ones(N,1);
rho = 10;
A = M'*M + rho*ones(p,p);
% Test whether A is invertible
if rcond(A) > SMALL % reciprocal condition number
IA = inv(A);
% ADMM parameters
lambda = 0.0002;
alpha = 1.1; % scaling factor
Iter = 30;
% Initialize abundance: sum of abundance equals to one
a_hat = rand(p, N);
for i = 1:p
a_hat(i, :) = a_hat(i, :) ./ sum(a_hat, 1);
end
z = a_hat;
u = zeros(p,N);
for iter = 1:Iter
x_tilde = z - u;
% Solve the QP problem with sum2one and positivity constraints
for i = 1:N
H = double(A);
f = double(-M'*y(:, i)-(1/rho)*x_tilde(:, i));
a_hat(:, i) = qpas(H, f, [], [], ones(p,p), ones(p,1), zeros(p,1), []); % C
end
z_tilde = a_hat + u;
% Denoising in 3D image domain
z_3d = permute(reshape(z_tilde,[p, sqrt(N), sqrt(N)]), [2,3,1]);
%% denoiser: bandwise
sigma = sqrt(lambda / rho);
z_3d_dn = zeros(size(z_3d)); % Initialize denoised 3D image
for i = 1:p
% z_3d_dn(:, :, i) = imnlmfilt(squeeze(z_3d(:, :, i))); % NLM
[~, z_3d_dn(:, :, i)]= BM3D(1, z_3d(:, :, i), sigma*255, 'lc'); % BM3D
% inPar.nSig=sqrt(lambda / rho);
% [z_3d_dn(:, :, i), outPar]= Denoiser(squeeze(z_3d(:, :, i)), ImgOrg, recMode{recMode_id}, inPar); %DnCNN
end
%% denoiser: 3D
% [z_3d_dn] = bm4d_1(1, z_3d, sigma); % BM4D
%% Converting into 2D domain
z = reshape(permute(z_3d_dn, [3, 1, 2]), [p, N]);
u = u + a_hat - z;
rho = alpha * rho;
end
end
%% Output
A_hat_graph = reshape(a_hat',100,100,p);
figure;
subplot(2, 2, 1);
imagesc(A_hat_graph(:,:,1),[0,1]);
axis off
subplot(2, 2, 2);
imagesc(A_hat_graph(:,:,2),[0,1]);
axis off
subplot(2, 2, 3)
imagesc(A_hat_graph(:,:,3),[0,1])
axis off
subplot(2, 2, 4)
imagesc(A_hat_graph(:,:,4),[0,1])
axis off