% demo1_SegSALSA

% Brief description
%
% This demo illustrates the use of SegSALSA algorithm introduced in
%
%   J. Bioucas-Dias, F. Condessa, and J. Kovacevic, "Alternating direction
%   optimization for image segmentation using hidden Markov measure field
%   models", in IS&T/SPIE Electronic Imaging, International Society for
%   Optics and Photonics,  pp. 90190P-90190P, February, 2014.
%
% with simulated data.
%
%%
% Summary:
%
%      The demo implements the following steps:
%
%      1) generate an image Y of labels by sampling  an MRF
%         Multilogistic Level prior
%
%      2) generate Gaussian feature vectors X(:,i), for i=1,...,n, given by
%
%                 X(:,i) = M(:,Y(i)) + N(:,i)
%
%         where M is a matrix holding the  class dependent mean   values and
%         N is Gaussian noise
%
%
%      3) For each pixel, compute the likelihood of each class.
%
%
%      4) Segment the  image using SegSALSA
%
%
% -----------------------------------------------------------------------
%
%
% Author: Jose Bioucas Dias, December, 2014

close all
clear all

rand('seed',314);
randn('seed',314);

%--------------------------------------------------------------------------
% generate an image of labels
%--------------------------------------------------------------------------
% % image dimension 
% dim=256;
% n = dim^2;
% % number of classes 
% no_classes = 3;
% % generate an MLL image
% beta=1.5;
% neigh = 2; % second order neighborhood
% generate classes
% Y = mll(dim,beta,no_classes,neigh);

%%
% -------------------------------------------------------------------------
% for faster execution, load precomputed image of classes
% -------------------------------------------------------------------------

no_classes = 3;
beta=1.5;
neigh = 2;
load classes_K3_n256_b1_5_neig2 Y

% ------------------------------------------------------------------------
% no_classes = 3;
% beta=3;
% neigh = 2; 
% load classes_K3_n256_b3_neig2 Y

% ------------------------------------------------------------------------
% no_classes = 10;
% beta=1.5;
% neigh = 2;
% load classes_K10_n256_b1_5_neig2 Y

% ------------------------------------------------------------------------
% no_classes = 10;
% beta=3;
% neigh = 2; 
% load classes_K10_n256_b3_neig2 Y
% ------------------------------------------------------------------------

[nr,nc] = size(Y);
n=nr*nc;

% confirm no_classes
no_classes = length(unique(Y(:)));

figure(10); figure(gcf)
imagesc(Y)
title('image of labels')
%%
% -------------- parameters for the feature vectors -----------
L = 1;     % dimension of the feature vectors

% standard deviation of the additive noise (component-wise)
sigma =  0.5;

                
%-------------------------------------------------------------------------
% generate data 
%--------------------------------------------------------------------------
M = 1: no_classes;
[L,p] = size(M);  % 


%--------------------------------------------------------------------------
% generate the data set
%--------------------------------------------------------------------------
Y = Y(:)';
X = zeros(L,n);
for i=1:no_classes
    np_class_i = sum(Y == i);
    X(:,Y == i) = M(i);
end


%--------------------------------------------------------------------------
% add Gaussian iid noise
%--------------------------------------------------------------------------
X = X + sigma*randn(size(X));

figure(30); figure(gcf)
imagesc(reshape(X,nr,nc))
title('Image of features')


%%
%--------------------------------------------------------------------------
% CLASSIFICATION
%--------------------------------------------------------------------------

% compute the likelihood images
P = zeros(no_classes,n);
for i=1:no_classes
    P(i,:) = 1/sqrt(2*pi)*exp(-(X-M(i)).^2/2/sigma^2);
end


[maxp,Y_ml] = max(P);
OA_ML =  sum(Y_ml==Y)/n

figure(50); figure(gcf)
imagesc(reshape(Y_ml,nr,nc));
st = sprintf('ML Classification: OA = %2.2f',OA_ML);
title(st)


%% segmentation using SegSALSA
 
 tau = 2.0;     % regularization parameter
 mu = 2.0;      % SALSA parameter 
 iters = 100;
 
 % segment the image
 figure(1000)
 P_segsalsa = imseg_hidd_vf(reshape(P',nr,nc,no_classes), tau,mu, iters);
 title('Latent probabilities for class 1 ')
 
 % remove boundary
 
[maxp,Y_segsalsa] = max(P_segsalsa);

Y_aux = reshape(Y_segsalsa, nr,nc);
Y_aux = Y_aux(:)';

OA_segsalsa =  sum(Y_aux==Y(:)')/n

figure(80); figure(gcf)
imagesc(reshape(Y_segsalsa,nr,nc));
st = sprintf('Segmentation: OA = %2.2f',OA_segsalsa);
title(st)




