function matches = link_neighbour_topological( re_part_A, re_part_B, R_search, K_feature, n_neighbours)
% matches = link_neighbour_topological(re_part_A, re_part_B, R_search, n_neighbors)
%
% Input:
%   re_part_A       3 x M   coordinates of particles in frame A
%   re_part_B       3 x N   coordinates of particles in frame B
%   R_search        1 x 1   radius of search for neighbouring particle patterns
%   K_feature		1 x 2   feature size in frames A and B
%						    i.e. number of neighbours to include in the search
%                           for local particle patterns
%	n_neighbours	1 x 1	number of neighbouring particles to search for
%
% Output:
%   matches         K x 2   indices of matching particle pairs (A->B)

%% sanity check
n_A         = size(re_part_A, 2);
n_B         = size(re_part_B, 2);
if nargin < 5
	n_neighbours = min(K_feature);
end
if n_A < K_feature(1)
    matches     = zeros(0, 2);
    return;
end
if n_B < K_feature(2)
    matches     = zeros(0, 2);
    return;
end
assert(size(re_part_A, 1) == 3);
assert(size(re_part_B, 1) == 3);
	
%% calculate feature arrays for fields A and B
% by default, use a feature of higher dimension in field A
% to account for the likely event that the features do not overlap
if numel(K_feature) == 1
	K_feature = [K_feature+1, K_feature];
else
	assert(K_feature(2) <= K_feature(1), 'feature dimension of frame B must be <= than A');
end
re_feature_A = build_feature(re_part_A, K_feature(1));
re_feature_B = build_feature(re_part_B, K_feature(2));

%% find the best match in B for particles in A
[inB, r_inB]= knnsearch(re_part_B.', re_part_A.', 'K', n_neighbours);
matches		= [1:n_A; nan(1,n_A)].';
G_best      = zeros(n_A,1);
no_match	= false(1,n_A);

for iA = 1 : n_A
	%% check whether nearest neighbours are within search distance
	ok			= r_inB(iA,:) < R_search;
	if ~any(ok)
		no_match(iA) = true;
		continue;
	end
	inB_ok		= inB(iA, ok);
	%% evaluate the likeness of candidate features in B to feature A
	G_candidates= feature_distance(re_feature_A(:,:,iA), re_feature_B(:,:, inB_ok));
	[G_best(iA),i_best]=min(G_candidates);
	matches(iA,2) = inB_ok(i_best);
end
% eliminate invalid matches
matches(no_match,:) = [];
G_best(no_match) = [];

%% ensure matches A->B are unique
% handle ties by picking the correspondence with the smallest feature
% distance
map				= sparse(matches(:,1), matches(:,2), 1./G_best, n_A, n_B);
[maxv,best]		= max(map,[],1);
iB				= find(maxv~=0);
iA				= best(maxv~=0); 

% assemble in array
matches			= [iA(:), iB(:)];

