function [feature_ids,scores]=cosmo_find_local_extrema(ds, nbrhood, varargin)
% find local extrema in a dataset using a neighborhood
%
% [feature_ids,scores]=cosmo_find_local_extrema(ds, nbrhood, ...)
%
% Inputs:
% ds dataset struct with one sample, i.e. size(ds.samples,1)
% must be 1. NaN values ds.samples are ignored.
% nbrhood neighborhood structure corresponding with dataset, with
% a field .neighbors so that .neighbors{k} contains the
% feature ids of neighbors of feature k.
% 'fitness', f optional function handle so that [score,id]=f(x)
% returns the index with the maximum 'feature score'.
% The default is @(x)max(x,[],2), i.e. a function that
% returns the maximum value and the index of the maximum
% value of the input.
% 'count', c optional number of elements to select. The default is
% Inf.
%
% Output:
% feature_ids 1xN vector with feature ids.
% scores 1xN vector with fitness scores. It holds for all I and
% J in 1:N that, if I<J, then:
% - fitness(ds.samples(I))>fitness(ds.samples(J))
% - intersect(nbrhood.neighbors{I},nbrhood.neighbors{J})
% is empty
% - if I_N=setdiff(nbrhood.neighbors{I},I), then
% fitness(ds.samples(I_N))<=fitness(ds.samples(I))
%
%
% Example:
% % Generate tiny dataset with 6 voxels
% ds=cosmo_synthetic_dataset('ntargets',1,'nchunks',1);
% cosmo_disp(ds.samples)
% %|| [ 2.03 0.584 -1.44 -0.518 1.19 -1.33 ]
% %
% % show unflattened shape
% cosmo_disp(squeeze(cosmo_unflatten(ds)))
% %|| [ 2.03 -0.518
% %|| 0.584 1.19
% %|| -1.44 -1.33 ]
% %
% nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
% % find local maxima within neighborhood of 1 voxel radius
% [feature_ids,scores]=cosmo_find_local_extrema(ds,nh);
% cosmo_disp(feature_ids)
% %|| [ 1 5 3 ]
% cosmo_disp(scores)
% %|| [ 2.03 1.19 -1.44 ]
% %
% % only return two feature ids
% [feature_ids,scores]=cosmo_find_local_extrema(ds,nh,'count',2);
% cosmo_disp(feature_ids)
% %|| [ 1 5 ]
% cosmo_disp(scores)
% %|| [ 2.03 1.19 ]
% %
% % use another fitness function, namely local minima
% [feature_ids,scores]=cosmo_find_local_extrema(ds,nh,'fitness',@min);
% cosmo_disp(feature_ids)
% %|| [ 3 4 ]
% cosmo_disp(scores)
% %|| [ -1.44 -0.518 ]
% %
% % find local maxima within neighborhood of 2 voxel radius
% nh=cosmo_spherical_neighborhood(ds,'radius',2,'progress',false);
% [feature_ids,scores]=cosmo_find_local_extrema(ds,nh);
% cosmo_disp(feature_ids)
% %|| [ 1 6 ]
% cosmo_disp(scores)
% %|| [ 2.03 -1.33 ]
%
%
% % A more sophisticated example showing how a table of local maxima
% % can be produced:
% %
% % Generate synthetic dataset with 700 voxels
% ds=cosmo_synthetic_dataset('size','big',...
% 'ntargets',1,'nchunks',1);
% %
% % ensure LPI orientation
% ds=cosmo_fmri_reorient(ds,'LPI');
% %
% % set minimum distance (in voxels) between local extrama (peaks)
% min_voxels_between_extrema=5;
% nh=cosmo_spherical_neighborhood(ds,...
% 'radius',min_voxels_between_extrema,...
% 'progress',false);
% %
% % find local maxima
% % (note: local minima are not found; this would require a
% % 'fitness' argument different from the default)
% [feature_ids,scores]=cosmo_find_local_extrema(ds,nh);
% %
% % define and apply threshold
% thr=1.96;
% msk=scores>thr;
% feature_ids=feature_ids(msk);
% scores=scores(msk);
% %
% % get coordinates of local extrama
% xyz=cosmo_vol_coordinates(ds,feature_ids);
% %
% % show table with coordinates and value at local extrema
% fprintf('x=% 3d y=% 3d z=% 3d score=%.3f\n',[xyz; scores]);
% %|| x= 11 y= 7 z= 3 score=3.679
% %|| x= -1 y= 3 z= -1 score=3.458
% %|| x= 27 y= 7 z= 7 score=3.361
% %|| x= 35 y= -1 z= 3 score=3.288
% %|| x= 23 y= 1 z= -1 score=2.899
% %|| x= 35 y= 11 z= -1 score=2.268
%
% Notes:
% - this function can be used to define regions of interest in a
% reproducible and objective manner.
% - to ignore particular features, set their value to NaN
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
