function [re_matched, im_matched, im_res, im_idx, c_idx] = match3d3c(cams, im_pos, opts)
	% [re_matched, im_matched, im_res, im_idx, c_idx] = match3d3c(cams, ij_pos, opts)
	%
	% given a set of particle locations observed in each of N cameras
	% perform stereo matching across pairs of cameras to find particles
	% that can be triangulated in 3 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
	%		gl_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_matched	3 x M array of positions of M matched particles
	% im_matched	2 x M x 3 array of image space coordinates of matched particles
	% im_dist		3 x M array of the residual distance on the image plane
	%				between particles and their back-projection on each of 3
	%				cameras
	% im_idx		3 x M array of the indices of particles in each camera
	% c_idx			3 x M array of indices of corresponding cameras
	
	%% basic parameters
	n_cameras			= length(cams);
	n_part				= cellfun(@(x) size(x,2), im_pos);

	%% convert to gpuArray
	%for c = 1 : n_cameras
	%	im_pos{c}		= gpuArray(im_pos{c});
	%end
	
	%% dewarp step
	im_pos_undistorted	= im_pos;
	for c = 1 : n_cameras
		im_pos_undistorted{c}= inverse_radial_distortion(im_pos{c}, cams(c).calib.ic, cams(c).calib.kr);
	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);
    ind_matches         = cell(1, n_stereo_pairs);
    n_pairs             = zeros(1, n_stereo_pairs);
	for k = 1 : n_stereo_pairs
		c1				= cpairs(k, 1);
		c2				= cpairs(k, 2);
		if ~opts.adjacency(c1,c2)
			% only compute matches for cameras which can see each other's
			% particles within the field of view
			continue;
		end
		
		[ic1,ic2]       = match_pairs(cams(c1), cams(c2), im_pos_undistorted{c1}, im_pos_undistorted{c2}, opts);
		ind_matches{k}  = [ic1(:), ic2(:)];
        n_pairs(k)      = numel(ic1);
        %fprintf('finding pairs C%d - C%d -> %d\n', c1, c2, numel(ic1));
    end
    
    %% find pairs that are consistent in three or more cameras
	% the result is stored in ij_idx, which is a M x N array
	% where M is the number of matches that are consistent amongst all three
	% cameras
    ctrips              = nchoosek(1:n_cameras, 3);
    n_camera_triplets   = size(ctrips, 1);
    re_matched_ary      = cell(1, n_camera_triplets);
    im_matched_ary      = cell(1, n_camera_triplets);
    im_res_ary          = cell(1, n_camera_triplets);
    im_idx_ary          = cell(1, n_camera_triplets);
    c_idx_ary           = cell(1, n_camera_triplets);
	dtype				= class(im_pos{1});
    
	for it = 1 : n_camera_triplets
        %% get combination of cameras
		c1              = ctrips(it, 1);
        c2              = ctrips(it, 2);
        c3              = ctrips(it, 3);
        if ~opts.adjacency(c1, c2) || ~opts.adjacency(c1, c3) || ~opts.adjacency(c1, c3)
            % only compute matches for cameras which share a field of view
            im_matched_ary{it}	= zeros(2,0,3,dtype);
			im_res_ary{it}		= zeros(3,0,dtype);
			re_matched_ary{it}  = zeros(3,0,dtype);
			im_idx_ary{it}		= zeros(3,0);
			c_idx{it}			= zeros(3,0);
			continue;
        end
        %% build adjacency matrix
        % find index of pair
        % this is guaranteed to work because c1 < c2 < c3 
        [~,k12]         = ismember([c1 c2], cpairs, 'rows'); 
        [~,k23]         = ismember([c2 c3], cpairs, 'rows');
        [~,k13]         = ismember([c1 c3], cpairs, 'rows');
		%A12				= sparse(ind_matches{k12}(:,1), ind_matches{k12}(:,2), 1, n_part(c1), n_part(c2));
		A23				= sparse(ind_matches{k23}(:,1), ind_matches{k23}(:,2), 1, n_part(c2), n_part(c3));
		A31				= sparse(ind_matches{k13}(:,2), ind_matches{k13}(:,1), 1, n_part(c3), n_part(c1)); % note transpose 
        %A11             = sparse(n_part(c1), n_part(c1));
        %A22             = sparse(n_part(c2), n_part(c2));
        %A33             = sparse(n_part(c3), n_part(c3));
        % construct block adjacency matrix
        %A123            = [A11 A12 A31.'; A12.' A22 A23; A31 A23.' A33];
        
        %% brute force the possible combinations that lead to closed cycles
		% adjacency matrices for pairs of cameras
		% start with all js accessible from is
		ic1                 = ind_matches{k12}(:,1);
        ic2                 = ind_matches{k12}(:,2);
        % now get potential paths from j->k, remembering which i we took
		n					= numel(ic2);
		[jn,kn,in]          = find(spdiags(ic1,0,n,n) * A23(ic2,:));
		jn					= ic2(jn);
		% finally, prune those that do not have a return loop
		has_loop			= A31(sub2ind([n_part(c3) n_part(c1)], kn, in)) > 0;
		in					= in(has_loop);
		jn					= jn(has_loop);
		kn					= kn(has_loop);
        in                  = in(:);
        jn                  = jn(:);
        kn                  = kn(:);
        
        %% triangulate and reject bad matches
        % triangulate
        im_matched			= cat(3, im_pos{c1}(:, in), im_pos{c2}(:, jn), im_pos{c3}(:, kn));
        [re_matched, ~, im_res] = triangulate(cams([c1 c2 c3]), im_matched, true);
        % identify bad matches
        im_res_ok           = max(im_res, [], 2, 'omitnan') < opts.R_tolerance_img;
        % reject bad matches
        lim_ok				= re_matched(1,:) >= opts.re_lim(1,1) ...
                            & re_matched(2,:) >= opts.re_lim(1,2) ...
                            & re_matched(3,:) >= opts.re_lim(1,3) ...
                            & re_matched(1,:) <= opts.re_lim(2,1) ...
                            & re_matched(2,:) <= opts.re_lim(2,2) ...
                            & re_matched(3,:) <= opts.re_lim(2,3);
        ok					= im_res_ok(:) & lim_ok(:);
        n_ok                = nnz(ok);
        
        %% accumulate in array
        % whilst rejecting bad matches
        re_matched_ary{it}  = re_matched(:, ok);
        im_matched_ary{it}  = im_matched(:, ok, :);
        im_res_ary{it}      = im_res(ok, :).';
        im_idx_ary{it}      = [in(ok,1) jn(ok,1) kn(ok,1)].';
        c_idx_ary{it}       = repmat([c1; c2; c3], 1, n_ok);
    end
    
    %% accumulate cell arrays
    re_matched      = cell2mat(re_matched_ary); % 3 x n_matched
    im_matched      = cell2mat(im_matched_ary); % 2 x n_matched x 3
    im_res          = cell2mat(im_res_ary);     % 3 x n_matched
    im_idx          = cell2mat(im_idx_ary);     % 3 x n_matched
    c_idx           = cell2mat(c_idx_ary);      % 3 x n_matched
end	
	