function [iA, iB, re_match, re_res, im_res] = match_pairs(camA, camB, im_A, im_B, opts)
	% [iA, iB, gl_match, gl_dist, ij_dist] = 
	%		match_pairs(camA, camB, ij_A, ij_B, opts)
	%
	% use stereoscopic matching match pairs of particle images observed in 
	% different cameras.
	%
	% input:
	%	camA					calibration for camera A
	%	camB					calibration for camera B
	%	im_A					location of particles in camera A in (dewarped) image coords
	%	im_B					location of particles in camera B in (dewarped) image coords
	%	opts					structure containing parameters for matching
	%		R_tolerance_img		threshold for maximum distance between re-projected
	%							particle image and its measured location on image plane

	%
	% return values:
	%	iA						index of image coordinate in camera A
	%	iB						index of image coordinate in camera B
	%	re_match				3D location of particle center, calculated by
	%							minimising least-squares distance from rays
	%	re_res					distance of particle to ray, in both cams
	%	im_res					distance of particle to reprojected image, in both
	%							cams

	iA         = [];
    iB         = [];
    re_match   = zeros(3,0);
    re_res     = zeros(2,0);
    im_res     = zeros(2,0);
    
    %% check im_A/im_B empty?
    if isempty(im_A) || isempty(im_B)
        return;
    end

    %% identify a common point in both fields of view
    % nominally at the centre of the measurement volume
    [b_adjacent, re_AB_lim] = check_adjacency(camA, camB, opts.re_lim);
    if ~b_adjacent
        return;
    end

    %% reduce set of particles for matching
    % particles can only be matched if their object space position is within
    % within re_AB_lim(1,:) and re_AB_lim(2,:). find the convex hull of
    % their projection to identify the set of particles which may be
    % triangulated within the measurement volume
    re_corners   = [re_AB_lim([1 1 1 1 2 2 2 2], 1).';
                    re_AB_lim([1 1 2 2 1 1 2 2], 2).';
                    re_AB_lim([1 2 1 2 1 2 1 2], 3).'];
    im_corners_A = pinhole_model(re_corners, camA.calib);
    im_corners_B = pinhole_model(re_corners, camB.calib);
    hull_A       = convhull(im_corners_A.');
    hull_B       = convhull(im_corners_B.');
    % check if particle lies within projection
    inA          = inpolygon(im_A(1,:), im_A(2,:), im_corners_A(1,hull_A), im_corners_A(2,hull_A));
    inB          = inpolygon(im_B(1,:), im_B(2,:), im_corners_B(1,hull_B), im_corners_B(2,hull_B));
    % cut down to the set of particles which are matchable
    im_A         = im_A(:,inA);
    im_B         = im_B(:,inB);
	n_A          = size(im_A, 2);
	n_B			 = size(im_B, 2);
    % retain index in original array for later
    inA          = find(inA);
    inB          = find(inB);
    inA          = inA(:);
    inB          = inB(:);
    
    if n_A == 0 || n_B == 0
        return;
    end
    
	%% calculate set of epipolar planes
	% imagine a line in object space
	% x_P = n0.s + x0
	% this can define a set of epipolar planes for cameras A and B
	% where each plane passes through x_A, x_B (the two camera positions)
	% and x_P, the point on the line
	%
	% given a point in image A or B, we can calculate the epipolar plane it
	% must lie on, and therefore its intersection with the line s
	%
	% if two particles have similar s, then they may correspond to one
	% another, since they lie on close epipolar planes. So the epipolar line
	% from P's image on A must be close to of P's image on B.
	%
	% We can use this to accelerate pair matching. We calculate the
	% projection onto the line 
	re_common       = mean(re_AB_lim).';
	[s_A,s_B,re_n0] = back_project_plane(camA.calib, camB.calib, re_common, im_A, im_B);
    
    %% find the tolerance for matching
    % rays reconstructed should match to within opts.R_tolerance_img
    % we need to find a corresponding tolerance between planes
    % obtain this from jacobian
    [im_common_A, im_jac_A] = pinhole_model(re_common, camA.calib);
	[im_common_B, im_jac_B] = pinhole_model(re_common, camB.calib);
    scale_A     = norm(im_jac_A * re_n0);
    scale_B     = norm(im_jac_B * re_n0);
    s_tolerance = max(2 * opts.R_tolerance_img ./ [scale_A scale_B]);
    
	%% search for pairs within given distance on line
	% s corresponds (approximately) to the projected distance between
	% a particle's back projection along a ray and the line x = x0 + n*s
	% 
	% we want to find pairs of rays that match reasonably closely,
	% so we're looking for pairs of s values below the maximum permissible
	% distance in object space. 
	
	% rangesearch returns cell array iAB(j)		j = 1 .. n_B
	% each cell contains indices i in range 1 .. n_A
	% corresponding to the elements of s_A which are within given distance
	% of s_B(j). 
	% here, we extract this information into vectors iA and iB,
	% which satisfy 
	%	|s_A(iA) - s_B(iB)| < opts.gl_max_dist
	
	ind_match			= rangesearch(gather(s_A(:)), gather(s_B(:)), s_tolerance);
	n_BtoA				= cellfun(@(x) size(x,2), ind_match);
	iB					= repelem(1 : n_B, n_BtoA).';
    iA					= cell2mat(ind_match(:).').';
	%[iA,iB]			= gridsearch(s_A, s_B, s_tolerance);
	
	%% back project pairs of particles
	yA					= im_A(:, iA);
	yB					= im_B(:, iB);
	cams				= [camA, camB];
	yAB					= cat(3, yA, yB);
	[re_match, re_res, im_res]	= triangulate(cams, yAB, false);
	
	%% reject invalid matches 
	% we check the criteria:
	%	are the image residuals below the required tolerance?
	%	is the particle triangulated within the volume bounds?
	is_im_res_ok		= all(im_res < opts.R_tolerance_img, 2);
	is_inside_bb		= re_match(1,:) >= re_AB_lim(1,1) ...
						& re_match(2,:) >= re_AB_lim(1,2) ...
						& re_match(3,:) >= re_AB_lim(1,3) ...
						& re_match(1,:) <= re_AB_lim(2,1) ...
						& re_match(2,:) <= re_AB_lim(2,2) ...
						& re_match(3,:) <= re_AB_lim(2,3);
	keep				= is_im_res_ok(:) & is_inside_bb(:);
	
	re_match			= re_match(:, keep);
	iA					= inA(iA(keep,1), 1);
	iB					= inB(iB(keep,1), 1);
	re_res				= re_res(keep, :);
	im_res				= im_res(keep, :);
end