function [ ind_dect, detect_RMUSIC, err_C ] = RMUSIC( B,th,Pyort,Mp,sel_mat )
% The Robust MUSIC algorithm for endmember detection
% ------------------Input ----
% B: library
% th: the error boud
% Pyort: orthogonal projector to the complement of R(Us)
% Mp: the number of endemebers one wants to select
% sel_mat: selected materials, which is the true endmember index set
%------------------output ----
% ind_dect: detected endmember index
% detect_RMUSIC: 1: sucess; 0: failure

N = length(sel_mat);
nx = 256;

Lib = B;
b2 = Pyort*B;
b1 = B - b2;

b1 = sqrt(repmat(sum(b1.^2),nx,1));
b2 = sqrt(repmat(sum(b2.^2),nx,1));

x = repmat(linspace(0,1,nx)',1,size(Lib,2));
% solve the optimization problem
[aux index] = max(  (b1 + th*x) ./ abs(b2 - th*(1-x).^0.5));
err_C = sqrt(1./ (1 + aux.^2));

[sort_err_C ind_errC]=sort(err_C,'ascend');
ind_errCC = sort(ind_errC(1:Mp),'ascend');
LibC = Lib(:,ind_errCC);                         % take out Mp endmembers with Mp smallest endmembers


% check whether we correctly detect the signatures
ind_dect = ind_errCC(1:Mp);
correct_set = sel_mat;
inter_set = intersect(ind_dect, correct_set);
inter_set = sort(inter_set,'ascend');

detect_RMUSIC=0;
if length(inter_set)==N
    if inter_set == sel_mat;
        detect_RMUSIC=1;
    end
else
    detect_RMUSIC=0;
end


end

