function [gl_pos, gl_dist, ij_dist, ij_idx, ij_pos] = match3d_scan(cams, im_pos_ary, opts)
	% gl_pos, gl_dist, ij_dist = match3d_scan(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 particle
	% pairs that match in all N cameras
	%
	% matching is augmented using depth information from scanning laser
	% sheet
	%
	%
	% inputs:
	% cams						Nx1 array of camera calibration structures
	%
	% ij_pos						Nx1 cell array of undistorted particle coordinates identified in each
	%								camera. each cell n contains a 3 x Np(n) array of
	%								coordinates of particles identified in that image
	%
	% opts						structure containing options for stereo matching
	%		scan_tolerance		threshold for maximum distance between particles
	%								in terms of number of sheets
	%		obj_tolerance		threshold for maximum distance between particle
	%								location and back-projected ray
	%		img_tolerance		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
	%
	% outputs:
	% gl_pos						3 x M array of positions of M matched particles
	% gl_dist					M x N array of the distance between particles and
	%								their back-projected rays in each camera 
	% ij_dist					M x N array of the distance on the image plane
	%								between particles and their back-projection
	% ij_idx						M x N array of the indices of each particle image
	%								location in ij_pos
	
	%% basic parameters
	n_cameras				= length(cams);
	n_part					= cellfun(@(x) size(x,2), im_pos_ary);
	
	if n_cameras > 2
		error('Scanning particle matching only implemented for 2 cameras');
	end
	camA						= cams(1);
	camB						= cams(2);
	
	%% get positions in dewarped coordinates
	im_pos_A					= inverse_radial_distortion(im_pos_ary{1}(1:2,:), camA.calib.ic, camA.calib.kr);
	im_pos_B					= inverse_radial_distortion(im_pos_ary{2}(1:2,:), camB.calib.ic, camB.calib.kr);
	sheet_pos_A				= im_pos_ary{1}(3,:);
	sheet_pos_B				= im_pos_ary{2}(3,:);
	
	int_pos_A				= round(sheet_pos_A);
	int_pos_B				= round(sheet_pos_B);
	sheet_min				= min([int_pos_A, int_pos_B]);
	sheet_max				= max([int_pos_A, int_pos_B]);
	n_sheets					= sheet_max - sheet_min + 1;
	sheet_range				= sheet_min : sheet_max;
	
	%% find matching particle pairs 
	% find pairs that match sheet #
	[iA_sheet, iB_sheet] = gridsearch(sheet_pos_A, sheet_pos_B, opts.scan_tolerance);
	adj_sheet				= sparse(iA_sheet, iB_sheet, 1, n_part(1), n_part(2));
	
	% match pairs for each sheet
	gl_pos_ary				= cell(1, n_sheets);
	gl_dist_ary				= cell(1, n_sheets);
	ij_dist_ary				= cell(1, n_sheets);
	ij_idx_ary				= cell(1, n_sheets);
	n_match					= zeros(1,n_sheets);
	for s = 1 : n_sheets
		%% find matches
		sheet_pos			= sheet_range(s);
		iA_cand				= find(int_pos_A == sheet_range(s));
		iB_cand				= find(ismember(int_pos_B, sheet_pos+(-2:2)));
		[ii, jj, gl_pos, gl_dist, ij_dist] = match_pairs(camA, camB, im_pos_A(:,iA_cand), im_pos_B(:,iB_cand), opts);
		iA						= iA_cand(1,ii);
		iB						= iB_cand(1,jj);
		
		%% reject those outside tolerance
		b_ok					= abs(sheet_pos_A(iA) - sheet_pos_B(iB)) < opts.scan_tolerance;
		iA						= iA(1,b_ok);
		iB						= iB(1,b_ok);
		gl_pos_ary{s}		= gl_pos(:, b_ok);
		gl_dist_ary{s}		= gl_dist(b_ok, :)';
		ij_dist_ary{s}		= ij_dist(b_ok, :)';
		ij_idx_ary{s}		= [iA; iB];
		n_match(s)			= nnz(b_ok);
	end
	
	%% collect arrays
	ok							= n_match > 0;
	gl_pos					= cell2mat(gl_pos_ary(ok));
	gl_dist					= cell2mat(gl_dist_ary(ok))';
	ij_dist					= cell2mat(ij_dist_ary(ok))';
	ij_idx					= cell2mat(ij_idx_ary(ok));
	
	% tack sheet position onto ij pos
	if any(n_match)
		ij_pos				= cat(3, im_pos_ary{1}(:, ij_idx(1,:)), ...
											im_pos_ary{2}(:, ij_idx(2,:)));
	else
		ij_pos				= zeros([3 0 2]);
	end
	
	%% reject non-unique matches
	if opts.b_unique
		n_used_A				= accumarray(ij_idx(1,:)', 1, [n_part(1) 1])';
		n_used_B				= accumarray(ij_idx(2,:)', 1, [n_part(2) 1])';
		ok_A					= find(n_used_A == 1);
		ok_B					= find(n_used_B == 1);
		ok						= ismember(ij_idx(1,:), ok_A) ...
								& ismember(ij_idx(2,:), ok_B);
		gl_pos				= gl_pos(:,ok);
		gl_dist				= gl_dist(ok,:);
		ij_dist				= ij_dist(ok,:);
		ij_idx				= ij_idx(:,ok);
		ij_pos				= ij_pos(:,ok,:);
	end
end
	