% Matlab demo code for the technical report with title  “Semisupervised 
% hyperspectral image unmixing using a variational Bayes algorithm”, 
% National Observatory of Athens, June, 2014
% 
% Copyright K. Themelis, A. Rontogiannis, K. Koutroumbas

clc ; clear; close all; format compact; format long eng;

load cuprite_ref.mat

Nen = 14; % number of endmembers
Np = size(x,2); %  number of pixels

% endmember matrix
Phi = x(:,randi(Np,1,Nen)); % select random pixels (bad)??? use VCA instead
Phi_gram = Phi.' * Phi;

[M,N] = size(Phi);

const = 1e-10;
MaxIter = 10000;

w_vb = zeros(N,Np);

for n = 1 : Np
    
    % VB estimation 
    w_vb(:,n) = mpLaplace_truncatedGaussian_vB(Phi, x(:,n), Phi_gram, MaxIter, const);
    
end

% reshape variables to display images
im = zeros(Lines, Columns, L); 
w_vb_im = zeros(Lines, Columns, Nen); 

for i = 1 : Lines, 
    for j = 1 : Columns, 
        im(i,j,:) = x(:,(j-1) * Lines + i ); 
        w_vb_im(i,j,:) = w_vb(:,(j-1) * Lines + i ); 
    end,
end

for i = 1 : Nen, figure; imagesc(w_vb_im(:,:,i),[0 1]); axis('image'); axis('off'); end
