function [re_part, D_part, P_part] = sipr_remove(re_part, D_part, P_part, tol)
	% [re_part, D_part] = sipr_merge(re_part, D_part, P_part tol)
	% merge particles which are too close together (below tol) 
	
	%% identify particles which have no nearest neighbour in search radius
	[~, d_nn]	= knnsearch(re_part.', re_part.', 'k', 2);
    if isempty(d_nn)
        return;
    end
	b_ok			= d_nn(:,2).' > 2*tol;
	
	%% isolate set which is OK
	re_ok			= re_part(:,b_ok);
	D_ok			= D_part(1,b_ok);
	P_ok			= P_part(1,b_ok);% a scalar property
	%% merge set which is not OK
	[re_merge, D_merge, P_merge] = remove(re_part(:,~b_ok), D_part(1,~b_ok), P_part(1,~b_ok), tol);
	%% concatenate
	re_part		= cat(2, re_ok, re_merge);
	D_part		= cat(2, D_ok,  D_merge);
	P_part		= cat(2, P_ok,  P_merge);
end

function [re_part, D_part, P_part] = remove(re_part, D_part, P_part, tol)
	% merge algorithm
	n_part					= size(re_part,2);
	if n_part < 1
		return;
	end
	
	%% find particles
	[iax,ibx,xsq]			= gridsearch(re_part(1,:), re_part(1,:), tol);
	[iay,iby,ysq]			= gridsearch(re_part(2,:), re_part(2,:), tol);
	[iaz,ibz,zsq]			= gridsearch(re_part(3,:), re_part(3,:), tol);
	
	%% create correspondence
	xsq						= sparse(iax, ibx, xsq, n_part, n_part);
	ysq						= sparse(iay, iby, ysq, n_part, n_part);
	zsq						= sparse(iaz, ibz, zsq, n_part, n_part);
	b_candidate				= sparse(iax, ibx, true(size(iax)), n_part, n_part) ...
								& sparse(iay, iby, true(size(iay)), n_part, n_part) ...
								& sparse(iaz, ibz, true(size(iaz)), n_part, n_part);
	rsq						= xsq + ysq + zsq;
	
	%% identify particles that are close together
	% the particle that is retained is brighter than all its neighbours
	b_close					= b_candidate & rsq < tol^2;
	[ii,jj]					= find(b_close);
	D2mat						= sparse(ii, jj, D_part(jj).^2, n_part, n_part);
	[~,i_best]				= max(D2mat, [], 2);
	b_best					= i_best == (1:n_part)'; % flag indicating if this particle is best amongst neighbours

	%% keep only the brightest
	re_part					= re_part(:,b_best);
	D_part					= D_part(1,b_best);
	P_part					= P_part(1,b_best);
end
	