function [pts] = remove_neighbourhood(pts, N)
IND = [];
for i=1:size(pts,1)
    d = sqrt(sum((pts - ones(size(pts,1),1)*pts(i,:)).^2,2));
    ind = find(d<N);
    ind = ind(find(ind>i));
    IND = [IND; ind];
end
pts(IND,:) = [];