%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% test particle pair matching with serial track
addpath(genpath([pwd '/../SerialTrack2']));
addpath(([pwd '/../../misc/mtimesx']));
addpath(([pwd '/../../mex/ndsyev']));

%% input parameters
delta_t				= 10e-3;
src_file			= 'L:\array\tcf\rescaled\z001_t001.mat';
n_neighbour_search	= 3;
n_subset			= 1e5;
algorithm			= 'myneightopo';

%% load particles
part		= load(src_file);

% evaluate subset
subset		= 1 : min(n_subset, part.n_part);
part.x		= part.x(:, subset);
part.u		= part.u(:, subset);
part.n_part = length(subset);

% evaluate displacement
re_part_t0	= part.x - part.u*delta_t;
re_part_t1	= part.x;
re_part_t2  = part.x + part.u*delta_t;

% randomise
P2			= randperm(part.n_part);
P0			= randperm(part.n_part);
re_part_t2	= re_part_t2(:, P2);
re_part_t0	= re_part_t0(:, P2);

% evaluate maximum displacement
u_norm		= sqrt(sum(part.u.^2, 1));
DeltaX_max	= max(u_norm)*delta_t*1.5;

%% try to match using serialtrack
% match forward step
fprintf('%d particles\n', part.n_part);
fprintf('matching frames 1 & 2...\n');
tic;
if strcmp(algorithm, 'serial')
	% SerialTrack algorithm
	m12			= f_track_serial_match3D(re_part_t1.', re_part_t2.', ...
						'f_o_s', DeltaX_max, ...
						'n_neighborsMax', n_neighbour_search);
elseif strcmp(algorithm, 'neightopo')
	m12			= f_track_neightopo_match3_fast(re_part_t1.', re_part_t2.', ...
						DeltaX_max, ...
						n_neighbour_search);
elseif strcmp(algorithm, 'myneightopo')
	m12			= link_neighbour_topological(re_part_t1, re_part_t2, ...
						DeltaX_max, ...
						[n_neighbour_search+2 n_neighbour_search]);
elseif strcmp(algorithm, 'nn')
	m12			= link_nearest_neighbour(re_part_t1, re_part_t2, DeltaX_max);
elseif strcmp(algorithm, 'hist')
	m12			= f_track_hist_match3(re_part_t1.', re_part_t2.', ...
						DeltaX_max, ...
						n_neighbour_search);
else
	error('unrecognised algorithm');
end
	
t_match		= toc;
fprintf('   took %0.3fs\n', t_match);

%% evaluate quality of matches
n_m12_error = nnz(m12(:,1) ~= P2(1,m12(:,2)).');
n_m12		= size(m12, 1);
fprintf('   found %d matches (%d errors)\n', n_m12, n_m12_error);

%%

