function [iX, iY, sqdist] = gridsearch(X, Y, R)
	% [iX, iY, sqdist] = gridsearch(X, Y, R)
	%
	% one dimensional search for neighbours within distance R
	%
	% find all elements in X which are within distance R of an element
	% in Y. returns iX, iY, which are 1xN arrays containing the indices of
	% matching pairs of elements in X and Y which are within distance R
	%
	
	NX				= numel(X);
	NY				= numel(Y);
	W				= R;
	X				= X(:);
	Y				= Y(:);
	
	if NX < 1 || NY < 1
		iX			= [];
		iY			= [];
		sqdist	= [];
		return;
	end
	
	%% Quantize elements in X and Y onto a grid of spacing R
	% if an element in X is within R of an element of Y, then it is
	% guaranteed to be adjacent or at the same grid index. However, it may
	% be up to 2R away, so we will have to reject such points later.
	
	Y_min			= min(Y);	
	Y_bin			= floor((Y - Y_min)/W) + 2;
	X_bin			= floor((X - Y_min)/W) + 2;
	
	%% create sparse map of bins within distance R of given Y index
	Y_bin_max	= max(Y_bin);
	iY				= (1 : NY)';
	map			= sparse([Y_bin-1; Y_bin; Y_bin+1], ...
								[iY     ; iY   ; iY     ], true, Y_bin_max+1, NY);
	
	%% match elements in X
	% indices of X values which may match those in Y (are not out of range)
	X_ok			= 1 <= X_bin & X_bin <= Y_bin_max+1;
	iX_ok			= find(X_ok);	
	% for those bins 
	[ii,iY]		= find(map(X_bin(X_ok), :));
	iX				= iX_ok(ii);
	
	%% reject X,Y pairs larger than R
	sqdist		= (X(iX(:)) - Y(iY(:)))'.^2;
	ok				= sqdist <= R^2;
	iX				= iX(ok);
	iY				= iY(ok);
	sqdist		= sqdist(ok);
end
