% // ====================================================================
% // This file is part of the Endmember Induction Algorithms Toolbox for MATLAB 
% // Copyright (C) Grupo de Inteligencia Computacional, Universidad del 
% // País Vasco (UPV/EHU), Spain, released under the terms of the GNU 
% // General Public License.
% //
% // Endmember Induction Algorithms Toolbox is free software: you can redistribute 
% // it and/or modify it under the terms of the GNU General Public License 
% // as published by the Free Software Foundation, either version 3 of the 
% // License, or (at your option) any later version.
% //
% // Endmember Induction Algorithms Toolbox is distributed in the hope that it will
% // be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
% // of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
% // General Public License for more details.
% //
% // You should have received a copy of the GNU General Public License
% // along with Endmember Induction Algorithms Toolbox. 
% // If not, see <http://www.gnu.org/licenses/>.
% // ====================================================================
%
%% [E,C] = EIA_CCA(data,p,t)
% 
% Manuel Grana <manuel.grana[AT]ehu.es>
% Miguel Angel Veganzones <miguelangel.veganzones[AT]ehu.es>
% Grupo de Inteligencia Computacional (GIC), Universidad del Pais Vasco /
% Euskal Herriko Unibertsitatea (UPV/EHU)
% http://www.ehu.es/computationalintelligence
% 
% Copyright (2011) Grupo de Inteligencia Computacional @ Universidad del Pais Vasco, Spain.
%
% Convex Cone Analysis (CCA) endmembers induction algorithm.
% ------------------------------------------------------------------------------
% Input:   data      : column data matrix [nvariables x nsamples]
%          p         : number of endmembers to be induced. If not provided it is calculated by HFC method with tol=10^(-5).
%          t         : tolerance for numerical errors. By default 10^(-6)
%
% Output:  E         : set of induced endmembers [nchannels x p]
%          C         : induced endmembers indexes vector [nsamples] with {0,1} values, where '1' indicates that the corresponding sample has been identified as an endmember. Some of the algorithms do not select pixels as the endmembers and, in that case C is empty.
%
% Bibliographical references:
% [1] Ifarraguerri, A., “Multispectral and hyperspectral image analysis with convex cones”, Geoscience and Remote Sensing, IEEE Transactions on, vol. 37, nº. 2, págs. 756-770, 1999.
function [E,C] = EIA_CCA(data,p,t)

%% Parameters
if (nargin < 1)
    error('Insufficient parameters');
end
if (nargin < 2 || p < 0)
    p = EIA_HFC(data,10^(-5));
end
if (nargin < 3 || t < 0)
    t = 10^(-6);
end

%% data size
[nvariables,nsamples] = size(data);

%% normalization
for i=1:nsamples
    norm_data = norm(data(:,i));
    data(:,i) = data(:,i) / norm_data;
end

%% Correlation matrix and eigen decomposition
correlation=data*data';
[V,D]=eig(correlation);
V=V(:,nvariables-p+1:nvariables);

%% Algorithm
% num of selected endmembers
endmembers_selec=0;
% init endmembers set
E = zeros(nvariables,p);
C = [];
% algorithm
while endmembers_selec<p
    % boolean vector indicating which bands are being selected
    selector=zeros(nvariables,1);
    % select p-1 bands randomly
    nselected=0;
    while nselected<p-1
        i=floor(rand*nvariables)+1;
        if ~selector(i)
            selector(i)=1;
            nselected=nselected+1;
        end
    end
    selector=find(selector);
    % equation system
    P=zeros(p-1);
    p1=zeros(p-1,1);
    for i=1:p-1
        P(i,:)=V(selector(i),1:p-1);
        p1(i)=V(selector(i),p);
    end
    A=P\p1;
    x = V*[A;1];
    if min(x)>-t
    	endmembers_selec=endmembers_selec+1
    	E(:,endmembers_selec)=x;
    end
end
