source: hypertools/hypertools/gldbm_td.m

%GLDBM_TD GLDB using top-down criterion (Fisher separability)
% 
%  W = GLDBM_TD(A)
% 
% INPUT  
%  A    Training dataset
%
% OUTPUT
%  W    best bases mapping
%
% DESCRIPTION 
% Generalised Best Bases mapping is a feature extraction technique for
% spectral datasets with two classes.  This function implements the
% top-down approach starting from a single group with all the wavelengths,
% recursively identifying the best possible split points.  Groups of
% adjacent features (wavelengths) are identified in the input dataset such
% that Fisher criterion is maximized in each group. Output dataset is
% created by projecting wavelengths in each group to the 1D Fisher
% subspace.  The procedure is defined for two class problems, if a
% multi-class dataset is supplied, a cell array of mappings will be
% computed.  Multi-class classifier, based on this routine, is implemented
% as PREXP algorithm ALG_GLDBM
%
% REFERENCE 
% Kumar,Ghosh,Crawford: Classification of Hyperspectral Data using
% Best-bases Feature Extraction Algorithm
%
% SEE ALSO
% DATASETS, MAPPINGS, GLDBM

% Copyright: P. Paclik, pavel@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: gldbm_td.m,v 1.5 2005/02/10 15:24:07 pavel Exp $

function [w,crit] = gldbm_td(a,arg2)

        % No arguments given: return untrained mapping.

        if (nargin < 1) | (isempty(a))
                w = mapping(mfilename);
                w = setname(w,'Best Bases mapping (top-down)');
                return
        end

        if (~exist('arg2','var'))     % Second argument is not a mapping: training.

                isdataset(a);           % Assert that A is a dataset.

                [m,k,c] = getsize(a); M = classsizes(a);
                
                if c>2
                        % for multi-class problem, create all pair-wise extractors
                        w={};
                        iter=[];
                        for i=1:c-1
                                for j=i+1:c
                                        fprintf(1,'(%d-%d',i,j);
                                        t=seldat(a,[i j]);
                                        eval(['w{i,j}=' mfilename '(t);']);
                                        fprintf(1,') ');
                                end
                        end
                        return
                end

                % counter with ten steps
                t=floor(k/10); % usually, only 80% bands is processed
                counts=t*ones(1,10);
                counts(1:mod(m,10))=counts(1:mod(m,10))+1;
                counts=cumsum(counts);
                
                Q=corrcoef(+a);
                [U,G]=meancov(a);

                m=0;              % level
                Nm=size(a,2);
                
                % in top-down approach, each group of bands is defined by its limits (l,u),
                % flag saying if it should be processed further and a criterion value:
                
                l=1;  % lower bounds of group-bands
                u=Nm; % upper bounds of group-bands
                process=1; % should the group be processed?
                I=0;    % criterion values for each group

                % compute crit for each single band group:
%               I=crit_eval(l,u,Q,+U(1,:),+U(2,:),G(:,:,1),G(:,:,2));

                while 1
                    
                    % we'll go over all the groups and process the processable ones:
                    temp=find(process==1);
                    gid=temp(1); % the 1st group in queue to be processed
%                   fprintf(1,' l=%d,u=%d ',l(gid),u(gid));
                    
                    % make possible merges:
                    % pl and pu refer to the possible left subgroup!
                    gc=u(gid)-l(gid)+1; % number of possible left subgroups
                    pl=repmat( l(gid), 1, gc-1 ); % possible beginnings of the left new subgroup
                    pu=pl(1):u(gid)-1;

                    % eval criterium for left subgroups:
                    crit1=crit_eval(pl,pu,Q,+U(1,:),+U(2,:),G(:,:,1),G(:,:,2));
                    
                    % eval criterion for right subgroups:
                    crit2=crit_eval(pu+1,repmat(u(gid),1,gc-1),Q,+U(1,:),+U(2,:),G(:,:,1),G(:,:,2));                
                    
                    maxval=max([crit1; crit2]);
                    [maxval,ind]=max(maxval);

                    newcrit=[crit1(ind) crit2(ind)];
                    
                    % pu(ind) is the best possible splitting index

                    % decide which of the two new subgroups should be processed further
                    newprocess=[0 0]; % by default: don't process
                    if crit1(ind)>I(gid) & ind>1
                        % decompose the left subgroup
                        newprocess(1)=1;
                    end
                    if crit2(ind)>I(gid) & gc-ind>1
                        % decompose the left subgroup
                        newprocess(2)=1;
                    end

                    % update high-level group definition:
                    if gid==1
                        % we're splitting the firrst group
                        l=[ 1       pu(ind)+1 l(2:end)];
                        u=[ pu(ind) u(1)      u(2:end)];
                        process=[newprocess process(2:end)];
                        I=[newcrit I(2:end)];
                    elseif gid==Nm
                        % we're splitting the last group
                        l=[ l(1:end-1) l(gid) pu(ind)+1];
                        u=[ u(1:end-1) pu(ind) Nm];
                        process=[process(1:end-1) newprocess];
                        I=[I(1:end-1) newcrit];
                    else
                        % we're in the middle
                        l=[ l(1:gid-1) l(gid) pu(ind)+1 l(gid+1:end)];
                        u=[ u(1:gid-1) pu(ind) u(gid)   u(gid+1:end)];
                        process=[process(1:gid-1) newprocess process(gid+1:end)];
                        I=[I(1:gid-1) newcrit I(gid+1:end)];                    
                    end
                    fprintf(1,'(%d) ', sum(process));
                    
                    if isempty(find(process==1))
                        % no candidates for progress left, return what
                        % we have at the moment.
                        [crit,clf]=crit_eval(l,u,Q,+U(1,:),+U(2,:),G(:,:,1),G(:,:,2));

                        % parameters
                        pars.clf=clf; pars.l=l; pars.u=u; pars.info_iter=m;
                        w = mapping(mfilename,'trained',pars,(1:length(l))',k,length(pars.l));
                        w = setname(w,'Best Bases mapping (top-down)');
                        
                        return
                    end
                end % while
        else % Second argument is a a mapping: testing.

                w = arg2;                       
                pars = getdata(w); % Unpack the mapping.
                [m,k] = getsize(a); 

                % create an output dataset
                out=[];
                for i=1:length(pars.clf)
                        % get the subset of bands...
                        t=+a(:,pars.l(i):pars.u(i));
                        % ... and project it to 1D Fisher subspace
                        out=[out t*pars.clf(i).w];
                end
        
                w = setdat(a,out);
                w=setfeatlab(w,(1:size(w,2))');
        end

return



function [I,clf]=crit_eval(l,u,Q,mua,mub,sigmaa,sigmab)

        I=zeros(1,length(l));

        clf=[];
        for i=1:length(l)
                
                C=min(min( Q( l(i):u(i), l(i):u(i) ) ));
                
                W=0.5*( sigmaa( l(i):u(i), l(i):u(i) )+sigmab( l(i):u(i), l(i):u(i) ) );
                t=( mua( l(i):u(i) ) - mub( l(i):u(i) ) )';
                B=t*t';
                
                w=pinv(W)*t;
                D=(w'*B*w)/(w'*W*w);
                
                I(i)=C*D;
                
                clf(i).w=w;
        end