function [row,col,max_local] = findLocalMaximum(val,radius)
    % Determine the local maximum of a given value
    %
    % Author :: Vincent Garcia
    % Date   :: 09/02/2007
    %
    % INPUT
    % =====
    % val    : the NxM matrix containing values
    % radius : the radius of the neighborhood
    %
    % OUTPUT
    % ======
    % row       : the row position of the local maxima
    % col       : the column position of the local maxima
    % max_local	: the NxM matrix containing values of val on unique local maximum
    %
    % EXAMPLE
    % =======
    % [l,c,m] = findLocalMaximum(img,radius);
    



    % FIND LOCAL MAXIMA BY DILATION (FAST) /!\ NON UNIQUE /!\
    % mask = fspecial('disk',radius)>0;
    % val2 = imdilate(val,mask);
    % index = val==val2;
    % [row,col] = find(index==1);
    % max_local = zeros(size(val));
    % max_local(index) = val(index);


    % FIND UNIQUE LOCAL MAXIMA USING FILTERING (FAST)
    mask  = fspecial('disk',radius)>0;
    nb    = sum(mask(:));
    highest          = ordfilt2(val, nb, mask);
    second_highest   = ordfilt2(val, nb-1, mask);
    index            = highest==val & highest~=second_highest;
    max_local        = zeros(size(val));
    max_local(index) = val(index);
    [row,col]        = find(index==1);


    % FIND UNIQUE LOCAL MAXIMA (FAST)
    % val_height  = size(val,1);
    % val_width   = size(val,2);
    % max_local   = zeros(val_height,val_width);
    % val_enlarge = zeros(val_height+2*radius,val_width+2*radius);
    % val_mask    = zeros(val_height+2*radius,val_width+2*radius);
    % val_enlarge( (1:val_height)+radius , (1:val_width)+radius ) = val;
    % val_mask(    (1:val_height)+radius , (1:val_width)+radius ) = 1;
    % mask  = fspecial('disk',radius)>0;
    % row = zeros(val_height*val_width,1);
    % col = zeros(val_height*val_width,1);
    % index = 0;
    % for l = 1:val_height
    %     for c = 1:val_width
    %         val_ref = val(l,c);
    %         neigh_val  = val_enlarge(l:l+2*radius,c:c+2*radius);
    %         neigh_mask = val_mask(   l:l+2*radius,c:c+2*radius).*mask;
    %         neigh_sort = sort(neigh_val(neigh_mask==1));
    %         if val_ref==neigh_sort(end) && val_ref>neigh_sort(end-1)
    %             index          = index+1;
    %             row(index,1)   = l;
    %             col(index,1)   = c;
    %             max_local(l,c) = val_ref;
    %         end
    %     end
    % end
    % row(index+1:end,:) = [];
    % col(index+1:end,:) = [];


end