function [re_pos, re_dist, im_dist, im_idx, im_matches] = match3d(cams, im_pos, opts)
	% re_pos, re_dist, im_dist = match3d(cams, im_pos, opts)
	%
	% given a set of particle locations observed in each of N cameras
	% perform stereo matching across pairs of cameras to find particle
	% pairs that match in all N cameras
	%
	% inputs:
	% cams		Nx1 array of camera calibration structures
	%
	% im_pos	Nx1 cell array of undistorted particle coordinates identified in each
	%			camera. each cell n contains a 2 x Np(n) array of
	%			coordinate pairs of particles identified in that image
	%
	% opts					structure containing options for stereo matching
	%
	%		R_tolerance_obj	threshold for maximum distance between particle
	%						location and back-projected ray
	%		R_tolerance_img	threshold for maximum distance between re-projected
	%						particle image and its measured location on image plane
	%		re_origin		coordinate of a common point in object space
	%						that is re-projected near the center of both
	%						images
	%		re_lim			limits of where particles can be reconstructed	
	%		b_dewarp		if true, input is interpreted to be specified in
	%						distorted image coordinates. output will be
	%						returned in distorted image coordinetes also.
	%
	% outputs:
	% re_pos	3 x M array of positions of M matched particles
	% re_dist	M x N array of the distance between particles and
	%			their back-projected rays in each camera 
	% im_dist	M x N array of the distance on the image plane
	%			between particles and their back-projection
	% im_idx	M x N array of the indices of each particle image
	%			location in im_pos
	
	%% basic parameters
	n_cameras			= length(cams);
	n_part				= cellfun(@(x) size(x,2), im_pos);
	
	%% dewarp step
	b_dewarp				= isfield(opts, 'b_dewarp') && opts.b_dewarp;
	im_pos_orig			= im_pos;
	if b_dewarp
		for c = 1 : n_cameras
			im_pos{c}	= inverse_radial_distortion(im_pos{c}, cams(c).calib.ic, cams(c).calib.kr);
		end
	end
	
	%% find matching particle pairs in pairs of cameras
	% also compute adjacency matrices for pairs in cameras
	cpairs				= nchoosek(1:n_cameras, 2);
	n_stereo_pairs		= nchoosek(n_cameras, 2);
	iA					= cell(1, n_stereo_pairs);
	iB					= cell(1, n_stereo_pairs);
	adj					= cell(1, n_stereo_pairs);
	for k = 1 : n_stereo_pairs
		cA				= cpairs(k, 1);
		cB				= cpairs(k, 2);
		[iA{k},iB{k}]	= match_pairs(cams(cA), cams(cB), im_pos{cA}, im_pos{cB}, opts);
		adj{k}			= sparse(iA{k}, iB{k}, 1, n_part(cA), n_part(cB));
	end
	n_matches			= cellfun(@(x) size(x, 2), iA);
	
	%% find pairs that are consistent in three or more cameras
	% the result is stored in im_idx, which is a M x N array
	% where M is the number of matches that are consistent amongst all three
	% cameras
	if n_cameras == 3
		% brute force the possible combinations that lead to closed cycles
		% adjacency matrices for pairs of cameras
		A12				= adj{1};
		A23				= adj{3};
		A31				= adj{2}.';

		% start with all js accessible from is
		[i,j]				= find(A12);
		% now get potential paths from j->k, remembering which i we took
		n					= numel(j);
		[jn,kn,in]		= find(spdiags(i,0,n,n) * A23(j,:));
		jn					= j(jn);
		% finally, prune those that do not have a return loop
		has_loop			= A31(sub2ind([n_part(3) n_part(1)], kn, in)) > 0;
		in					= in(has_loop);
		jn					= jn(has_loop);
		kn					= kn(has_loop);
		% store
		im_idx			= [in, jn, kn];
	else
		error('Algorithm not implemented for four or more cameras');
	end
	
	%% get ij coordinates of consistent n-tuplets
	n_triplets			= size(im_idx, 1);
	im_matches			= zeros(2, n_triplets, n_cameras);
	for c = 1 : n_cameras
		im_matches(:,:,c) = im_pos_orig{c}(:, im_idx(:,c));
	end
	
	%% triangulate and reject bad matches
	% triangulate
	[re_pos, re_dist, im_dist] = triangulate(cams, im_matches, b_dewarp);
	% identify bad matches
	im_ok				= max(im_dist, [], 2) < opts.R_tolerance_img;
	% check if outside reconstruction limits
	lim_ok				= re_pos(1,:) >= opts.re_lim(1,1) ...
						& re_pos(2,:) >= opts.re_lim(1,2) ...
						& re_pos(3,:) >= opts.re_lim(1,3) ...
						& re_pos(1,:) <= opts.re_lim(2,1) ...
						& re_pos(2,:) <= opts.re_lim(2,2) ...
						& re_pos(3,:) <= opts.re_lim(2,3);
	ok					= im_ok & lim_ok(:);
	% reject bad matches
	im_matches			= im_matches(:, ok, :);
	re_pos				= re_pos(:, ok);
	re_dist				= re_dist(ok, :);
	im_dist				= im_dist(ok, :);	
	im_idx				= im_idx(ok, :);
end	
	