function [tracks,tracks2f,re_unmatched] = track_frame3(camera, frame, mopts)
	% [tracks, tracks2f, re_unmatched] = track_frame3(frame, mopts)
	%	
	% track particles across three frames using time-symmetric neighbour tracking
	% algorithm
	%
	% Input:
	%   camera  1xNC array of camera model structures
	%
	%	frame	1x3 cell array, with structures containing the fields
	%		n_part		1x1			number of particles in this frame
	%		re_part		3xn_part	particle coordinates
	%
	%	mopts	options structure, containing the following fields
	%		R_search_obj	search radius for particle pairs, mm
	%		K_feature		dimension of features to track with NN-algorithm
	%		n_search		maximum number of particles to consider within
	%						search radius
	%		t				time coordinate 
	%		v_pred			predictor field function, v_pred(x) (optional)
	%
	% Output:
	%	tracks	1 x n_tracks  track array structure containing
	%		x		3x1		fitted position at t=0
	%		u		3x1		fitted velocity at t=0
	%		re_pos	3x3		measured object space particle position at t=1..3
	%		im_pos	2x3xNC	measured image space particle position at t=1..3
	%		im_pred	2x3xNC	predicted image position 
	%
	%	tracks2f		track array structure containing partially matched
	%
	%	re_unmatched 1 x 3 cell array containing arrays
	%		3 x n_unmatched	positions of unmatched particles at t=1..3
	
	n_cameras	= numel(camera);
	dtype		= class(frame{1}.re_part);
	
	%% evaluate predictor field
	re_part_pred_t0	= cell(1, 3);
	if isfield(mopts, 'v_pred') && isa(mopts.v_pred, 'function_handle')
		% use predictor field to map particles back to t = 0
		for it = 1 : 3
			re_part_pred_t0{it} = frame{it}.re_part - mopts.t(it)*mopts.v_pred(frame{it}.re_part);
		end
		r_search = mopts.R_search_predictor;
	else
		% no predictor field used
		for it = 1 : 3
			re_part_pred_t0{it} = frame{it}.re_part;
		end
		r_search = mopts.R_search_obj;
	end
	
	%% frame 2->3
	fprintf('frame 2->3, matching %d->%d objects ... ', frame{2}.n_part, frame{3}.n_part);
	match23 = link_neighbour_topological(re_part_pred_t0{2}, ...
				re_part_pred_t0{3}, ...
				r_search, ...
				mopts.K_feature, ...
				mopts.n_search);
	n_match_23 = size(match23, 1);
	fprintf('%d correspondences\n', n_match_23);
	
	%% frame 2->1
	fprintf('frame 2->1, matching %d->%d objects ... ', frame{2}.n_part, frame{1}.n_part);
	match21 = link_neighbour_topological(re_part_pred_t0{2}, ...
				re_part_pred_t0{1}, ...
				r_search, ...
				mopts.K_feature, ...
				mopts.n_search);
	n_match_21 = size(match21, 1);
	fprintf('%d correspondences\n', n_match_21);
	
	%% frame 1->3
	if nargout >= 2
		fprintf('frame 1->3, matching %d->%d objects ... ', frame{1}.n_part, frame{3}.n_part);
		match13 = link_neighbour_topological(re_part_pred_t0{1}, ...
					re_part_pred_t0{3}, ...
					r_search, ...
					mopts.K_feature, ...
					mopts.n_search);
		n_match_13 = size(match13, 1);
		fprintf('%d correspondences\n', n_match_13);
	end
	
	%% identify set of particles where all 3 frames are linked
	[b_both,locb]= ismember(match21(:,1), match23(:,1));
	n_match_123 = nnz(b_both);
	match123	= [match21(b_both,[2 1]) match23(locb(b_both), 2)];
	%match123	= [match21(b_both,[2 1]) match23(locb(b_both), :)];
	%assert(all(match123(:,2) == match123(:,3)))
	fprintf('%d matched in both frames\n', n_match_123);
	
	%% evaluate particle velocities
	% positions at times 1, 2, & 3
	re_match_t1		= frame{1}.re_part(:, match123(:,1));
	re_match_t2		= frame{2}.re_part(:, match123(:,2));
	re_match_t3		= frame{3}.re_part(:, match123(:,3));
	[x0,u]			= track_fit_linear3(mopts.t, re_match_t1, re_match_t2, re_match_t3);
	
	%% check the image space residuals
	b_large_error	= zeros(3, n_match_123);
	im_pred_ary     = cell(1,3);
	im_meas_ary   	= cell(1,3);
	
	for it = 1 : 3
		%% project particles on all cameras at this timestep
		R			= frame{it};
		re_pred		= x0 + mopts.t(it)*u;
		im_pred_full= zeros(2, n_match_123, n_cameras, dtype);
		for ic = 1 : n_cameras
			im_pred_full(:,:,ic) = pinhole_model(re_pred, camera(ic).calib, true);
		end
		
		ip			= match123(:,it);
		im_meas_full= im_pred_full;
		for ic = 1 : 3
			js					= sub2ind([n_match_123 n_cameras], 1:n_match_123, R.c_idx(ic,ip));
			im_meas_full(:,js)	= R.im_matched(:,ip,ic);
		end
		im_res		= im_meas_full - im_pred_full;
		im_res_norm = sqrt(sum(im_res.^2,1));
		b_large_error(it,:) = any(im_res_norm > mopts.R_tolerance_img, 3);
		%% store
		im_meas_ary{it} = im_meas_full;
		im_pred_ary{it} = im_pred_full;
	end
	% eliminate tracks with large residuals from linear fitting procedure 
	b_reject		= any(b_large_error, 1);
	fprintf('%d rejected due to velocity discrepancy\n', nnz(b_reject));
	
	%% build track objects
	tracks			= cell(1, n_match_123);
	for ip = 1 : n_match_123
		t			= struct();
		t.u			= u(:,ip);
		t.x			= x0(:,ip);
		t.im_pred	= cat(2, im_pred_ary{1}(:,ip,:), im_pred_ary{2}(:,ip,:), im_pred_ary{3}(:,ip,:));
		t.im_pos	= cat(2, im_meas_ary{1}(:,ip,:), im_meas_ary{2}(:,ip,:), im_meas_ary{3}(:,ip,:));
		t.re_pos	= cat(2, re_match_t1(:,ip), re_match_t2(:,ip), re_match_t3(:,ip));
		tracks{ip}	= t;
	end
	tracks(b_reject)= [];
	tracks			= cell2mat(tracks);
	
	%% optional: build second set of tracks based on two frames only
	if nargout < 2
		return;
	end
	b_23			= ~ismember(match23(:,1), match123(~b_reject,2));
	b_21			= ~ismember(match21(:,1), match123(~b_reject,2));
	match23_ok		= match23(b_23,:);
	match21_ok		= match21(b_21,:);
	% matches from 2->3 only that we haven't already considered
	x0_23			= frame{2}.re_part(:, match23_ok(:,1));
	x1_23			= frame{3}.re_part(:, match23_ok(:,2));
	u_23			= (x1_23-x0_23)/(mopts.t(3)-mopts.t(2));
	% matches from 1->2 only that we haven't already considered
	x0_21			= frame{2}.re_part(:, match21_ok(:,1));
	x1_21			= frame{1}.re_part(:, match21_ok(:,2));
	u_21			= (x1_21-x0_21)/(mopts.t(1)-mopts.t(2));
	% matches from 1->3, none of these are discounted because there is no
	% particle to discount at time step 2
	x0_13			= frame{1}.re_part(:, match13(:,1));
	x1_13			= frame{3}.re_part(:, match13(:,2));
	u_13			= (x1_13-x0_13)/(mopts.t(3)-mopts.t(1));
	x_13			= (x0_13+x1_13)/2; % position at t=0 assumes t(1)=-t(3)
	assert(mopts.t(1) == -mopts.t(3));
	
	n_tracks_23		= nnz(b_23);
	n_tracks_21		= nnz(b_21);
	n_tracks_13		= n_match_13;
	tracks23		= cell(1, n_tracks_23);
	tracks21		= cell(1, n_tracks_21);
	tracks13		= cell(1, n_tracks_13);
	for ip = 1 : n_tracks_23
		tracks23{ip}= struct('x', x0_23(:,ip), 'u', u_23(:,ip));
	end
	for ip = 1 : n_tracks_21
		tracks21{ip}= struct('x', x0_21(:,ip), 'u', u_21(:,ip));
	end
	for ip = 1 : n_tracks_13
		tracks13{ip}= struct('x', x_13(:,ip), 'u', u_13(:,ip));
	end
	tracks2f		= cat(2, tracks23, tracks21, tracks13);
	tracks2f		= cell2mat(tracks2f);
	n_2f			= size(tracks2f,2);
	fprintf('%d particles in at least 2 frames\n', n_2f);

	%% project tracks which appear in two frames only
	x2f				= [tracks2f.x];
	u2f				= [tracks2f.u];
	im_pred_2f		= zeros(2, n_2f, n_cameras, 3, dtype);
	for it = 1 : 3
		re_pred		= x2f + mopts.t(it)*u2f;
		for ic = 1 : n_cameras
			im_pred_2f(:,:,ic,it) = pinhole_model(re_pred, camera(ic).calib, true);
		end
	end
	im_pred_2f		= permute(im_pred_2f, [1 4 3 2]);
	im_pred_2f		= mat2cell(im_pred_2f, 2, 3, n_cameras, ones(1, n_2f));
	[tracks2f.im_pred]=im_pred_2f{:};
	
	%% identify all unmatched particles
	if nargout < 3
		return;
	end
	re_unmatched	= cell(1,3);
	i_matched		= cell(1,3);
	i_matched{1}	= [match123(:,1); match13(:,1); match21(:,2)];
	i_matched{2}	= [match123(:,2); match23(:,1); match21(:,1)];
	i_matched{3}	= [match123(:,3); match13(:,2); match23(:,2)];
	for it = 1 : 3
		re_unmatched{it} = frame{it}.re_part;
		re_unmatched{it}(:,i_matched{it}) = [];
	end
end