% FUNCTION [ LOCAL_MAX_LIST, LOCAL_MAX_VALUE_LIST ] = LOCAL_MAX( DATA, NBEST )
%
% For each row of DATA, return the list of local maxima ( "greater or equal" meaning ).
%
% 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_FAST, LOCAL_MAX_STRICT, LOCAL_MAX_STRICT_FAST
%
% By Guillaume Lathoud 2004 - lathoud@idiap.ch

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

% Deal with parameters
%
% This is paranoid, in order to detect bugs in programs calling "local_max"
  
  if nargin < 1
    error( 'local_max: needs one input argument !' );
  end

  if nargout < 1
    error( 'local_max: needs one output argument !' );
  end
  
  if ~isnumeric( data )
    error( 'local_max: input argument "data" must be NUMERIC !' );
  end
  
  if any( imag( data(:) ) )
    error( 'local_max: input argument "data" must be REAL !' );
  end

  siz = size( data );
  
  if ~exist( 'nbest', 'var' )
    nbest = [];

  elseif ~isreal( nbest )
    error( 'local_max: input argument "nbest" must be a real, positive value' );
  
  elseif isnan( nbest )
    error( 'local_max: input argument "nbest" must be a real, positive value' );
  
  elseif nbest < 0
    error( 'local_max: input argument "nbest" must be a real, positive value' );
  
  elseif 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
    
  % 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 );
  
  if nargout > 1
    local_max_value_list = cell( prod( siz( 2:end ) ), 1 );
  end
  
  for a = 1:length( local_max_list )
    
    d = diff( [ -Inf; data( :, a ); -Inf ] );
    
    % Two ways, depending on "nbest" parameter
    
    if ~isempty( nbest )

      % Limit the number of local maxima
      
      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;
      
      if nargout > 1
	local_max_value_list{ a } = v;
      end
    
    else
    
      % Return all local maxima
      
      local_max_list{ a } = find( ( d( 1:end-1 ) >= 0 ) & ( d(2:end) <= 0 ) );
      if nargout > 1
	local_max_value_list{ a } = data( local_max_list{ a }, a );
      end
    
    end
    
  end
  
  % Reshape the result to match with "data" size
  
  local_max_list = reshape( local_max_list, [ siz( 2:end ) 1 ] );

  if nargout > 1
    local_max_value_list = reshape( local_max_value_list, [ siz( 2:end ) 1 ] );
  end