% FUNCTION [ LOCAL_MAX_LIST, LOCAL_MAX_VALUE_LIST ] = LOCAL_MAX_FAST( DATA, NBEST )
%
% For each row of DATA, return the list of local maxima ( "greater or
% equal" meaning ). This is the same code as "LOCAL_MAX", but most
% safety tests have been removed. Use it only when you are sure of
% your own code.
%
% DATA: a matrix of real values of any dimensionality D [ N1 N2 N3 ... ND ].
%
% LOCAL_MAX_LIST: a cell array of dimension [ N2 N3 ... ND ]. 
%   LOCAL_MAX_LIST{ I2, I3, ... , ID } is a vector.
%   This vector is the list of row indices of local maxima of DATA( :, I2, I3, ..., ID ). 
%
% LOCAL_MAX_VALUE_LIST: a cell array of dimension [ N2 N3 ... ND ]. 
%   LOCAL_MAX_VALUE_LIST{ I2, I3, ... , ID } is a vector.
%   This vector is the list of values of local maxima of DATA( :, I2, I3, ..., ID ). 
%   ( same order as in LOCAL_MAX_LIST )
%
% (optionally) NBEST: integer, a maximum bound on the number of local maxima 
%   returned in each vector.
%
% Note: if NBEST is not used, then the results are ordered by
% increasing row index. If NBEST is used, then the results are ordered
% by decreasing value (global maximum comes first, second best comes
% second, etc.).
%
% Therefore, if you want all local maxima, ordered by decreasing
% value, use NBEST = +Inf.
%
%
%
% See also: LOCAL_MAX, LOCAL_MAX_STRICT, LOCAL_MAX_STRICT_FAST
%
% By Guillaume Lathoud 2004 - lathoud@idiap.ch

function [ local_max_list, local_max_value_list ] = local_max_fast( data, nbest )

  siz = size( data );
  
  if nargout < 2
    
    % We work in flat 2-dimensional space,
    % and at the end we reshape "local_max_list".
    %
    % This way any dimensionality of "data" can be processed.
    
    local_max_list       = cell( prod( siz( 2:end ) ), 1 );
    
    % Two ways, depending on "nbest" parameter
    
    if nargin > 1

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Limit the number of local maxima

      if nbest == 0
	if ~strcmp( warning, 'off' )
	  disp( 'local_max: WARNING! nbest = 0' );
	end
	
	local_max_list = cell( siz( 2:end ) );
	
	if nargout > 1
	  local_max_value_list = cell( siz( 2:end ) );
	end
	
	return;
	
      end

      
      for a = 1:length( local_max_list )
	
	d = diff( [ -Inf; data( :, a ); -Inf ] );
	
	ind = find( ( d( 1:end-1 ) >= 0 ) & ( d(2:end) <= 0 ) );
	
	v = data( ind, a );
	[ v, ind2 ] = sort( -v );
	
	if length( ind ) > nbest
	  ind2 = ind2( 1:nbest );
	end
	
	ind = ind( ind2 );
	
	local_max_list{ a }       = ind;
	
      end
      
      % Reshape the result to match with "data" size
      
      local_max_list       = reshape( local_max_list, [ siz( 2:end ) 1 ] );
      
    else
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
      % Return all local maxima
      
      for a = 1:length( local_max_list )
	
	d = diff( [ -Inf; data( :, a ); -Inf ] );
	
	local_max_list{ a }       = find( ( d( 1:end-1 ) >= 0 ) & ( d(2:end) <= 0 ) );
	
      end      

      % Reshape the result to match with "data" size
      
      local_max_list        = reshape( local_max_list, [ siz( 2:end ) 1 ] );
      
    end

  else   % if nargout > 1
    
    % ----------------------------------------------------------------------
    % Also return values at the local maxima
    
    % We work in flat 2-dimensional space,
    % and at the end we reshape "local_max_list".
    %
    % This way any dimensionality of "data" can be processed.
    
    local_max_list       = cell( prod( siz( 2:end ) ), 1 );
    local_max_value_list = cell( prod( siz( 2:end ) ), 1 );
    
    % Two ways, depending on "nbest" parameter
    
    if nargin > 1

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Limit the number of local maxima

      if nbest == 0
	if ~strcmp( warning, 'off' )
	  disp( 'local_max: WARNING! nbest = 0' );
	end
	
	local_max_list = cell( siz( 2:end ) );
	
	if nargout > 1
	  local_max_value_list = cell( siz( 2:end ) );
	end
	
	return;
	
      end

      
      for a = 1:length( local_max_list )
	
	d = diff( [ -Inf; data( :, a ); -Inf ] );
	
	ind = find( ( d( 1:end-1 ) >= 0 ) & ( d(2:end) <= 0 ) );
	
	v = data( ind, a );
	[ v, ind2 ] = sort( -v );
	v = -v;
	
	if length( v ) > nbest
	  v   = v( 1:nbest );
	  ind2 = ind2( 1:nbest );
	end
	
	ind = ind( ind2 );
	
	local_max_list{ a }       = ind;
	local_max_value_list{ a } = v;
	
      end
      
      % Reshape the result to match with "data" size
      
      local_max_list       = reshape( local_max_list, [ siz( 2:end ) 1 ] );
      local_max_value_list = reshape( local_max_value_list, [ siz( 2:end ) 1 ] );
      
    else
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
      % Return all local maxima
      
      for a = 1:length( local_max_list )
	
	d = diff( [ -Inf; data( :, a ); -Inf ] );
	
	local_max_list{ a }       = find( ( d( 1:end-1 ) >= 0 ) & ( d(2:end) <= 0 ) );
	local_max_value_list{ a } = data( local_max_list{ a }, a );
	
      end      

      % Reshape the result to match with "data" size
      
      local_max_list        = reshape( local_max_list, [ siz( 2:end ) 1 ] );
      local_max_valuse_list = reshape( local_max_value_list, [ siz( 2:end ) 1 ] );
      
    end
    
  end