function reconnect(obj)
	% reconnect track objects using matching in position-velocity space
	% 
	% options for track reconnection are in recon structure
	% it has the following fields:
	%
	% 	n_search		search interval between ends of tracks, in frames
	% 	n_truncate		number of data points to truncate before reconnecting
	%	C_threshold		1 x n_tau vector of cost threshold values C(tau)
	%	tau				1 x n_tau time lag for each value in C_threshold
	% 	r_collision		minimum distance between particles
	%					particles closer than this distance are marked as
	%					ghosts
	%	f_cutoff		cutoff frequency for filtering tracks with
	%					smoothing spline
	%	
	
	%% parameters for track reconnection
	recon			= obj.opts.recon;
	% search interval
	recon.T_search	= diff(obj.t_axis([1 2])) * recon.n_search;

	%% fit smoothing spline to tracks
	obj.fit_spline(recon.f_cutoff);

	%% find track end points
	n_tracks= numel(obj.tracks);
	i_start = arrayfun(@(t) find(~t.ext,1,'first'), obj.tracks);
	i_end   = arrayfun(@(t) find(~t.ext,1,'last' ), obj.tracks);
	% get velocity at end points
	% (x1,u1,t1) is feature at end point
	% (x2,u2,t2) is feature at start point
	x1		= cell2mat(arrayfun(@(t,ind) t.re_pred(:,1,ind), obj.tracks, i_end  , 'UniformOutput', false));
	u1		= cell2mat(arrayfun(@(t,ind) t.re_vel(:,1,ind) , obj.tracks, i_end  , 'UniformOutput', false));
	t1		= obj.t_axis(1,i_end  );
	
	x2		= cell2mat(arrayfun(@(t,ind) t.re_pred(:,1,ind), obj.tracks, i_start, 'UniformOutput', false));
	u2		= cell2mat(arrayfun(@(t,ind) t.re_vel(:,1,ind) , obj.tracks, i_start, 'UniformOutput', false));
	t2		= obj.t_axis(1,i_start);	
	
	%% link tracks 
	% matches(:,1) is index of earlier track
	% matches(:,2) is index of later track
	matches  = link_position_velocity(x1,u1,t1,x2,u2,t2,recon);
	n_matches= size(matches,1);
	fprintf('reconnecting %d tracks\n', n_tracks);
	fprintf('found %d matches\n', n_matches);

	%% combine tracks together
	% merge tracks A and B so that A->AB, B->AB
	% we need to do this in chronological order
	% e.g.
	%	merging 1&2 then 2&3 results in
	% 		1&2->12, 2&3->123, 1->12, 2->123, 3->123 
	%	whereas 2&3 then 1&2 results 
	%		2&3->23, 1&2->123, 1->123, 2->123, 3->23 
	% we can avoid redundancy by eliminating the first track in each merge

	% sort matches chronologically
	t2_match	= t2(1, matches(:,2));
	reorder		= sort(t2_match, 'ascend');
	matches		= matches(reorder,:);
	
	% combine tracks
	fields		= fieldnames(obj.tracks);
	for n = 1 : n_matches
		% merge properties
		trk1	= obj.tracks(matches(n,1));
		trk2	= obj.tracks(matches(n,2));
		% merge flags
		trk1	= trk1.ext & trk2.ext;
		trk1	= trk1.oob | trk2.oob;
		% merge numeric properties
		if regexp(fields{n}, 're_.+|im_.+', 'once')
			prop1	= trk1.(fields{n});
			prop2	= trk2.(fields{n});
			prop1(:,:,~trk2.ext) = prop2(:,:,~trk2.ext);
			trk1.(fields{n}) = prop1;
		end
		% assign both tracks
		obj.tracks(matches(n,:)) = trk1;
	end

	% eliminate redundant tracks
	obj.tracks(matches(:,1)) = []; 

	%% fit smoothing spline to new tracks
	% interpolate values 
	obj.fit_spline(recon.f_cutoff);
	n_tracks = numel(obj.tracks);
	for n = 1 : n_tracks
		trk		= obj.tracks(n);
		i_start = find(~trk.ext, 1, 'first');
		i_end   = find(~trk.ext, 1, 'last');
		% interpolate position values where there is no data
		trk.re_pos(:,:,trk.ext) = trk.re_pred(:,:,trk.ext);
		% invalid values beyond track
		trk.re_pos(:,:,1:i_start-1) = nan;
		trk.re_pos(:,:,i_end+1:end) = nan;
	end
end

function D = pairdist(X1, X2)
	% D = pairdist(X1, X2)
	%
	% X1 is a 1 x 7 vector encoding:
	% position	x1 = X1(1:3)'
	% velocity  u1 = X1(4:6)'
	% time      t1 = X1(7)'
	%
	% X2 is a N x 7 matrix encoding the same for a set of points (x2,u2,t2)
	% we compute the generalised distance
	% 
	% D = |dx|^2 + |du|^2*|dt|^2/12
	% where
	%	dx = x2 - x1
	%	du = u2 - u1
	%	dt = t2 - t1

	x1 = X1(1,1:3);
	u1 = X1(1,4:6);
	t1 = X1(7);
	x2 = X2(:,1:3);
	u2 = X2(:,4:6);
	t2 = X2(:,7);

	dx	 = x2-x1;
	dudt = (u2-u1).*(t2-t1);
	D    = sum(dx.^2 + dudt.^2, 2) / 12;
end

