function [particles, matches] = match(obj, prt_ary)
    % [particles, matches] = match(obj, prt_ary)
    %
    % 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_c cameras
	%
	% Input:
    %   prt_ary     1 x n_c cell array
    %       contains structure arrays of particle centers found in images
    %       with the fields
    %           centroid    2 x 1   image space centroid
    %           length      2 x 1   (optional) length/orientation of rod
    %           radius      1 x 1   (optional) radius of sphere
    %           intensity   1 x 1   (optional) intensity of particle
    %
    % Output:
    %   particles   1 x n_match structure array with the following fields
    %       re_pos      3 x 1   object space position
    %       re_length   3 x 1   object space length
    %       re_radius   1 x 1   object space radius
    %       im_pos      2 x n_c measured image space position
    %       im_length   2 x n_c measured image space length
    %       im_radius   2 x n_c measured image space radius
    %       im_intens   1 x n_c measured image space intensity
    %       im_proj     2 x n_c projected image space position
    %       im_res      1 x n_c residual in image space position
    %       im_idx      1 x n_c index of particle in prt_ary
    %
    %   matches     n_match x n_c array with indices of 2d particles in
    %       prt_ary used to reconstruct 3d particles
    
    %% basic parameters
    camera              = obj.setup.camera;
	n_cameras			= length(camera);
	n_part				= cellfun(@(x) size(x,2), prt_ary);
	topts				= obj.opts.triangulation;

	%% extract particle properties into separate arrays
	im_centroid_ary		= cell(1, n_cameras);
	im_length_ary		= cell(1, n_cameras);
	im_radius_ary		= cell(1, n_cameras);
	im_intens_ary		= cell(1, n_cameras);
    b_has_length        = false;
    b_has_radius        = false;

	for c = 1 : n_cameras
		im_centroid_ary{c}	= [prt_ary{c}.centroid];
		% additional properties for rods/spheres
		if isfield(prt_ary{c}, 'length')
			im_length_ary{c}    = [prt_ary{c}.length];
            b_has_length        = true;
		else
			im_length_ary{c}	= nan(2, n_part(c));
		end
		if isfield(prt_ary{c}, 'radius')
			im_radius_ary{c}    = [prt_ary{c}.radius];
            b_has_radius        = true;
		else
			im_radius_ary{c}    = nan(1, n_part(c));
		end
		if isfield(prt_ary{c}, 'intensity')
			im_intens_ary{c}    = [prt_ary{c}.intensity];
		else
			im_intens_ary{c}	= nan(1, n_part(c));
		end
	end

	%% matching on particle centers
	topts.re_lim		= obj.setup.gl_bb_lim;
	[re_pos,~,~,matches]= match3d(camera, im_centroid_ary, topts);
	n_match				= size(re_pos, 2);

	% find reprojection
	im_pos				= zeros(2, n_match, n_cameras);
    im_proj             = zeros(2, n_match, n_cameras);
	im_length			= zeros(2, n_match, n_cameras);
	im_radius			= zeros(1, n_match, n_cameras);
	im_intens			= zeros(1, n_match, n_cameras);
	im_jac              = zeros(2, 3, n_match, n_cameras);
    for c = 1 : n_cameras
        % measured position
        im_pos(:,:,c)   = im_centroid_ary{c}(:, matches(:,c));
        % reprojected position and jacobian
        [im_proj(:,:,c), im_jac(:,:,:,c)] = pinhole_model(re_pos, camera(c).calib);
        % length
        im_length(:,:,c)= im_length_ary{c}(:, matches(:,c));
        % radius
        im_radius(:,:,c)= im_radius_ary{c}(1, matches(:,c));
        % intensity
        im_intens(:,:,c)= im_intens_ary{c}(1, matches(:,c));
    end
	% image space residual, 1 x n_match x n_cameras
	im_res				= sqrt(sum((im_proj-im_pos).^2,1));

	%% find properties of rod-like particles
    % for a rod of size dx, the image space projection dy_i is related by
    % the jacobian J_i as
    %   J_i * dx = dy_i 
    % we set up a system of equations
    %   [J_1; J_2; J_3; ...] * dx = [dy_1; dy_2; dy_3; ...]
    % and solve for the unknown length dx
    
    if b_has_length
        % get the length of rods in image space as a 2C x 1 x n_match array
        % and the jacobian as a 2C x 3 x n_match array
        im_length_perm  = reshape(permute(im_length, [1 3 2]), [2*n_cameras 1 n_match]);
        im_jac_perm     = reshape(permute(im_jac,  [1 4 2 3]), [2*n_cameras 3 n_match]);

        % the image space measurement of length is ambiguous in its
        % orientation. therefore, consider every possible orientation of
        % particles in image space, and pick the combination which has the
        % lowest residual
        n_perms         = 2^(n_cameras-1);
        power           = (1:n_cameras-1)';
        power           = repelem(power,2);
        sign_perm       = mod(0:n_perms-1, 2.^power) >= 2.^(power-1);
		sign_perm		= [ones(2,n_perms); sign_perm];
        sign_perm       = sign_perm*2-1;
        
        % evaluate residual for every combination
        re_length_perm  = zeros(3, 1, n_match, n_perms);
        im_length_resnorm= zeros(n_match, n_perms);
        for p = 1 : n_perms 
			% solve in least squares sense
            y						= sign_perm(:,p).*im_length_perm;
            re_length_perm(:,:,:,p) = pagemldivide(im_jac_perm, y);
            im_length_res			= pagemtimes(im_jac_perm, re_length_perm(:,:,:,p)) - y;
            im_length_resnorm(:,p)	= sum(im_length_res.^2, 1);
        end
        % pick the optimal combination
        [~,i_best]      = min(im_length_resnorm, [], 2);
		k_best			= sub2ind([n_match n_perms], 1 : n_match, i_best.');
		re_length		= re_length_perm(:,k_best);
        %re_length       = squeeze(re_length(:,:,:,i_best)); % 3 x n_match
    else
        re_length       = nan(3, n_match);
    end
    
    %% find image space properties of spherical particles
    % for a sphere of radius R, the image space projection r_i is 
    % related to the jacobian as
    %    r_i = maximum of R * |J_i * e| for |e| = 1
    % i.e. 
	%    r_i = R * sigma_i 
	% where sigma_i is the largest singular value of J_i
	if b_has_radius
        im_JtJ			= pagemtimes(im_jac, 'transpose', im_jac, 'none');
        sigma_m			= sqrt(ndsyev(single(im_JtJ), 's-'));
        sigma_m			= double(sigma_m(1,:,:));

        % solve this system r = X*R in the least squares sense for R
        % where X = [sigma_1; sigma_2; ...]
        sigma_m			= permute(sigma_m, [3 1 2]);
        im_radius_perm  = permute(im_radius, [3 1 2]);

        XtX				= pagemtimes(sigma_m, 'transpose', sigma_m, 'none');
        Xtr				= dot(im_radius_perm, sigma_m, 1);
        re_radius		= Xtr ./ XtX;
        re_radius		= reshape(re_radius, [1 n_match]);
    else
        re_radius       = nan(1, n_match);
    end
    
	%% collect properties into arrays
    flatten         = @(x) reshape(x, [1 numel(x)]);
    num2cell_flat   = @(x, dim) flatten(num2cell(x,dim));
	particles		= ...
		struct(	're_pos',		num2cell_flat(re_pos   , 1), ...
				're_length',    num2cell_flat(re_length, 1), ...
				're_radius',    num2cell_flat(re_radius, 1), ...
                'im_idx',       num2cell_flat(matches, 2), ...
				'im_pos',		num2cell_flat(permute(im_pos   , [1 3 2]), [1 2]), ...
				'im_proj',	    num2cell_flat(permute(im_proj  , [1 3 2]), [1 2]), ...
				'im_res',		num2cell_flat(permute(im_res   , [1 3 2]), [1 2]), ...
				'im_length',	num2cell_flat(permute(im_length, [1 3 2]), [1 2]), ...
				'im_intens',	num2cell_flat(permute(im_intens, [1 3 2]), [1 2]), ...
				'im_radius',	num2cell_flat(permute(im_radius, [1 3 2]), [1 2]));

	%   particles   1 x n_match structure array with the following fields
    %       re_pos      3 x 1   object space position
    %       re_length   3 x 1   object space length
    %       re_radius   1 x 1   object space radius
    %       im_pos      2 x n_c measured image space position
    %       im_length   2 x n_c measured image space length
    %       im_radius   2 x n_c measured image space radius
    %       im_intens   1 x n_c measured image space intensity
    %       im_proj     2 x n_c projected image space position
    %       im_res      1 x n_c residual in image space position
end
	
