function  pairs= pairing_mesh(surf,value,L,sigma)


% function  pairs= pairing_mesh(surf, value, L, sigma)
%
% Construct the min-max pairing of function value along a sphere.
%
% surf  : spherical mesh
% value : functiond defined on surface mesh coordinates
% L     : degree of weighted-SPHARM expansion
% sigma : bandwidth
%
%
% (C) Moo K. Chung 2009-
%
% mkchung@wisc.edu
%
% Department of Statistics and Medical Informatics
% University of Wisconsin-Madison
%
% Update history: Feb 2, 2009 initial code
%               July 24, 2010 works for toy examples
%



% find all local min and max.

tri=surf.faces;
coord=surf.vertices;

[nbr, degree] = FINDnbr(tri);

n_vertex=max(max(tri));
n_tri=size(tri,1);

maxind=[];
minind=[];

% THIS CODE WORKS BUT PRODUCES TOO MANY CRITICAL POINTS THAT ARE LIKELY TO BE NOISE
%for i=1:n_vertex
%    flag=sum(value(i) > value(nbr(i,1:degree(i))));
    % check if the node value is greater than neighbors
    % flag = node degree only if it is local maximum
%    if flag== degree(i)
 %       maxind=[maxind i];  % node index for all local max
%    end;
    % checking local minimum is done similarly
%    flag=sum(value(i) < value(nbr(i,1:degree(i))));
%    if flag == degree(i)
%        minind=[minind i];  % node index for all local min
%    end
%end;



% smooth data by putting value into the 1st coordinate


coord=zeros(size(surf.vertices,1),3);
coord(:,1)=value;

input.vertices=coord;
input.faces=surf.faces;

[output,fourier]= SPHARMsmooth2(input,surf, L,sigma);

square=SPHARM2square(fourier.x(:,:,1),L,sigma);
[lmax, lmin] = figure_extrema(square);

% using maximums and minimums, we reconstruct triangulation
% this simplifies the construction of persistence diagram substantially

%coordmin=coord(minind,:);
%coordmax=coord(maxind,:);
%surfmin.vertices=coordmin;
%surfmax.vertices=coordmax;

%[thetamin,varphimin]=EULERangles(surfmin);  
%[thetamax,varphimax]=EULERangles(surfmax);   

thetamin=lmin.I; varphimin=lmin.J;
thetamax=lmax.I; varphimax=lmax.J;
 
% hold on;
% plot(100*varphimin,100*thetamin,'+k');
% hold on
% plot(100*varphimax,100*thetamax,'or');

x=[thetamin ;thetamax];
y=[varphimin; varphimax];


trim= delaunay(x,y);
%hold on; triplot(trim,y,x ,'--k');
%print -r300 -dtiff persistence3

[nbrm, degm] = FINDnbr(trim);

% min max indicator 
%minmax=[-ones(length(thetamin),1); ones(length(thetamax),1)];
%indm=[minind'; maxind'];
%critical =value(indm);
%g =value(minind); %local minimums
%h= value(maxind); %local maximums

% all critical values
g=lmin.value;
h=lmax.value;
critical = [g , h];

gsort=sort(g);
hsort=sort(h);
hind=length(g)+1:length(g) + length(h);

m=length(g);
pair=[];


for i=m:-1:1
    ind_c=find(critical==gsort(i));
    nbr_ind=nbrm(ind_c,1:degm(ind_c));
    nbr_ind=intersect(nbr_ind, hind);
    
    if ~isempty(nbr_ind)
        nbr_value =critical(nbr_ind);
        argmin=find(nbr_value>gsort(i));
        if ~isempty(argmin)
            hstar = min(nbr_value(argmin));

            hstar_ind=find(critical==hstar);
            hind=setdiff(hind,hstar_ind);
            pair=[pair; [gsort(i) hstar]];
        end;
%     else
%         % case where the smallest max is taken by other min
%         
%         argmin=find(critical(hind)>gsort(i));
%         hstar = min(critical(argmin))
% 
%         hstar_ind=find(critical==hstar);
%         hind=setdiff(hind,hstar_ind);

    end;
    
end;

pairs=pair;
return

%plot(pair(:,1),pair(:,2),'.')

