function OUT = neighbors(L,k)

%NEIGHBORS computes the values of the kth nearest neighbors to points.

%WRITTEN BY: Tony E. Smith, 2/22/01

%INPUTS: L = (Nx2) matrix of (Xi,Yi) = coordinates
%        
%   		k = order of the nearest neighbors to be considered.

%        
%NOTE: The program produces only one kth neighbor (the first one found)
%      and hence is not appropriate for regular grid data.
%
%OUTPUTS: OUT = (Nx4)matrix specified.
%         OUT(i,1)   = row number of kth nearest neighbor to i
%         OUT(i,2:3) = coordinates of kth nearest neighbor to i
%			 OUT(i,4)   = kth nearest neighbor distance to i


N = length(L(:,1)) ;

if max(k,N)==k 
   
   error('k is too big for data set');    
     
end


p = k + 3 ;

for i = 1:N
   
   D = sqrt( (L(i,1)-L(:,1)).^2 + (L(i,2) - L(:,2)).^2 ); %distances to i
   
   D(i)= max(D) + 1 ; %Eliminate diagonal elements as n-neighbors.   
    
   Y(:,1) = [1:N]' ;
   
   Y(:,2) = D ;
   
   Y = sortrows(Y,2) ; %orders distances to i
   
   nn_k = Y(k,1);   
   
   OUT(i,1) = nn_k;
   
   OUT(i,2:3) = L(nn_k,:);
   
   OUT(i,4) = D(nn_k);
   
     
end