function [ D, Di, Q, beta ] = DRE( D, tol, ind )
% DRE Dictionary Redundancy Elimination
%
% With help of column pivoting it tries  to choose  a stable basis from a given
% dictionary.  
%  
% Usage:    [ Dnew, Di, Q, beta ] = DRE( D, tol, ind ); 
%           [ Dnew, Di ] = DRE( D );
%
% Inputs:
%   D       dictionary of normalized atoms 
%   tol     tolerance set for considering as linearly dependent atom [default tol= 
%           1.0000e-7]
%   ind     (optional) indices determining the initial subspace 
%   []      can be used for  tol and ind
%
% Outputs:
%   D       sub-dictionary extracted from D containing linearly independent 
%           atoms
%   Di      indices of atoms in new D w.r.t to original D 
%   Q       orthonormal functions spanning the same space as new D 
%   beta    biorthogonal functions to new D
%
% References: 
%   L. Rebollo-Neira, "Dictionary redundancy elimination", IEE Proceedings - Vision, Image
%   and Signal Processing, Vol(151,1), 31-34 (2004).
%
% See also DREOr DREOrp Biorthog.

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/

name='DRE';  %name of routine
% setting of default tolerance parameter
if (nargin<3) |(isempty(ind==1)); ind=[];  end
if (nargin<2) |(isempty(tol==1)); tol=1e-7;end
%original number of function in  dictionary D
N=size(D,2); 
Q=D;
beta=[];
Di=1:N;     % index vector of full-dictionary
numind=length(ind);

%===============================
% Main algorithm: at kth iteration
%===============================
for k=1:N
  c=zeros(1,N);
  %column pivoting for chosen atoms
  if k<=numind 
    for p=1:length(ind),
      [test,z]=ismember(ind(p),Di);
      if test~=1 error('Demanded index %d is out of dictionary',ind(p));end
      c(p)=norm(Q(:,z));
    end
    [max_c,q1]=max(c);
    if max_c<tol 
      error('Demanded atom with index %d is dependent on the others',ind(q1));
      break; 
    end
    [test,q]=ismember(ind(q1),Di);ind(q1)=[];
  else 
    %column pivoting
    for p=k:N, c(p)=norm(Q(:,p)); end
    [max_c,q]=max(c);
    if max_c<tol  k=k-1;break; end %stopping criterion
  end  
  if q~=k
    Q(:,[k q])=Q(:,[q k]); % swap k-th and q-th columns
    D(:,[k q])=D(:,[q k]);
    Di([k q])=Di([q k]);
  end	
  %re-orthogonalization of Q(:,k) w.r.t Q(:,1),..., Q(:,k-1)
  if k>1
    for p=1:k-1
      Q(:,k)=Q(:,k)-(Q(:,p)'*Q(:,k))*Q(:,p);
    end
   end
  nork=norm(Q(:,k)); 
  Q(:,k)=Q(:,k)/nork; %normalization
  temp=Q(:,k)/nork; %this will be beta(:,k)
  %another re-orthogonalization 
  for i=1:k-1
    temp=temp-Q(:,i)*(Q(:,i)'*temp);
  end   
  %compute biorthogonal functions beta from 1 to k-1
  if k>1
   beta=beta-temp*(D(:,k)'*beta);
  end
  %kth biorthogonal function
  beta(:,k)=temp;
  %orthogonalization of Q(:,k+1:n) w.r.t Q(:,k)
  if k==N, break; end   
  for p=k+1:N
    % orthogonalization of Q(:,p) wrt Q(:,k)
    Q(:,p)=Q(:,p)-(Q(:,k)'*Q(:,p))*Q(:,k);
  end
end
D=D(:,1:k); %new sub-dictionary
Q=Q(:,1:k); %orthonormal functions
Di=Di(1:k); % array of indices
fprintf('\n** %s - Dictionary Redundancy Elimination **\n',name);
fprintf('From %d atoms in this dictionary %d atoms are linearly independent.\n',N,k);
fprintf('The tolerance parameter was set up tol=%g.\n',tol); 

% error testing
% orthogonality
ErrorTest(Q);
% biorthogonality
ErrorTest(D,beta);

%Last Modification Laura Rebollo-Neira 2009.
%
%Copyright (C) 2006 Miroslav ANDRLE and Laura REBOLLO-NEIRA
%
%This program 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 2 of the License, or (at your option) any later version.
%
%This program 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 this program;
%if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
%Boston, MA  02110-1301, USA.
