Code covered by the BSD License

# Simple Tracker

### Jean-Yves Tinevez (view profile)

01 Dec 2011 (Updated )

Simple multiple particle tracker with gap closing.

simpletracker(points, varargin)
```function [ tracks adjacency_tracks A ] = simpletracker(points, varargin)
% SIMPLETRACKER  a simple particle tracking algorithm that can deal with gaps
%
% *Tracking* , or particle linking, consist in re-building the trajectories
% of one or several particles as they move along time. Their position is
% reported at each frame, but their identiy is yet unknown: we do not know
% what particle in one frame corresponding to a particle in the previous
% frame. Tracking algorithms aim at providing a solution for this problem.
%
% |simpletracker.m| is - as the name says - a simple implementation of a
% tracking algorithm, that can deal with gaps. A gap happens when one
% particle that was detected in one frame is not detected in the subsequent
% one. If not dealt with, this generates a track break, or a gap, in the
% frame where the particule disappear, and a false new track in the frame
% where it re-appear.
%
% first created between each frame pair, using by default the hungarian
% found to be the closest (euclidean distance). By virtue of the hungarian
% algorithm, it is ensured that the sum of the pair distances is minimized
% over all particules between two frames.
%
% Then a second iteration is done through the data, investigating track
% ends. If a track beginning is found close to a track end in a subsequent
% track, a link spanning multiple frame can be created, bridging the gap
% and restoring the track. The gap-closing step uses the nearest neighbor
%
% INPUT SYNTAX
%
% tracks = SIMPLETRACKER(points) rebuilds the tracks generated by the
% particle whose coordinates are in |points|. |points| must be a cell
% array, with one cell per frame considered. Each cell then contains the
% coordinates of the particles found in that frame in the shape of a
% |n_points x n_dim| double array, where |n_points| is the number of points
% in that frame (that can vary a lot from one frame to another) and |n_dim|
% is the dimensionality of the problem (1 for 1D, 2 for 2D, 3 for 3D,
% etc...).
%
% tracks = SIMPLETRACKER(points, KEY, VALUE, ...)  allows to pass extra
% parameters to configure the tracking process settings. Accepted KEYS &
% VALUES are:
%
% 'Method' - a string, by default 'Hungarian'
% Specifies the method to use for frame-to-frame particle linking. By
% default, the hungarian method is used, which ensures that a global
% optimum is found for each frame pair. The complexity of this algorithm is
% in O(n^3), which can be prohibitive for problems with a large number of
% point in each frame (e.g. more than 1000). Therefore, it is possible to
% use the nearest-neighbor algorithm by setting the method to
% 'NearestNeighbor', which only achieves a local optimum for a pair of
% points, but runs in O(n^2).
%
% 'MaxLinkingDistance' - a positive number, by default Inifity.
% Defines a maximal distance for particle linking. Two particles will not
% be linked (even if they are the remaining closest pair) if their distance
% is larger than this value. By default, it is infinite, not preventing nay
%
% 'MaxGapClosing' - a positive integer, by default 3
% Defines a maximal frame distance in gap-closing. Frames further way than
% this value will not be investigated for gap closing. By default, it has
% the value of 3.
%
% 'Debug' - boolean flag, false by default
% Adds some printed information about the tracking process if set to true.
%
% OUTPUT SYNTAX
%
% track = SIMPLETRACKER(...) return a cell array, with one cell per found
% track. Each track is made of a |n_frames x 1| integer array, containing
% the index of the particle belonging to that track in the corresponding
% frame. NaN values report that for this track at this frame, a particle
% could not be found (gap).
%
% Example output: |track{1} = [ 1 2 1 NaN 4 ]| means that the first track
% is made of the particle 1 in the first frame, the particule 2 in the
% second frame, the particle 1 in the 3rd frame, no particle in the 4th
% frame, and the 4th particle in the 5th frame.
%
% [ tracks adjacency_tracks ] = SIMPLETRACKER(...) return also a cell array
% with one cell per track, but the indices in each track are the global
% indices of the concatenated points array, that can be obtained by
% |all_points = vertcat( points{:} );|. It is very useful for plotting
% applications.
%
% [ tracks adjacency_tracks A ] = SIMPLETRACKER(...) return the sparse
% between a source particle (row) and a target particle (column) where
% there is a 1. Rows and columns indices are for points in the concatenated
% points array. Only forward links are reported (from a frame to a frame
% later), so this matrix has no non-zero elements in the bottom left
% diagonal half. Reconstructing a crude trajectory using this matrix can be
% as simple as calling |gplot( A, vertcat( points{:} ) )|
%
% VERSION HISTORY
%
% * v1.0 - November 2011 - Initial release.
% * v1.1 - May 2012 - Solve memory problems for large number of points.
%                   - Considerable speed improvement using properly the
%                   sparse matrices.
%                   - Use the key/value pair syntax to configure the
%                   function.
% * v1.3 - August 2012 - Fix a severe bug thanks to Dave Cade
%
% Jean-Yves Tinevez < jeanyves.tinevez@gmail.com> November 2011 - 2012

%% Parse arguments

p = inputParser;
defaultDebug                = false;
defaultMaxGapClosing        = 3;
defaultMethod               = 'Hungarian';
expectedMethods = { defaultMethod, 'NearestNeighbor' };

@(x) any(validatestring(x, expectedMethods)));

p.parse( varargin{:} );

debug                   = p.Results.Debug;
max_gap_closing         = p.Results.MaxGapClosing;
method                  = p.Results.Method;

if debug
fprintf('Frame to frame linking using %s method.\n', method);
end

n_slices = numel(points);

current_slice_index = 0;
row_indices = cell(n_slices, 1);
column_indices = cell(n_slices, 1);
unmatched_targets = cell(n_slices, 1);
unmatched_sources = cell(n_slices, 1);
n_cells = cellfun(@(x) size(x, 1), points);

if debug
fprintf('%03d/%03d', 0, n_slices-1);
end

for i = 1 : n_slices-1

if debug
fprintf(repmat('\b', 1, 7));
fprintf('%03d/%03d', i, n_slices-1);
end

source = points{i};
target = points{i+1};

switch lower(method)

case 'hungarian'
[target_indices , ~, unmatched_targets{i+1} ] = ...

case 'nearestneighbor'
[target_indices , ~, unmatched_targets{i+1} ] = ...

end

unmatched_sources{i} = find( target_indices == -1 );

% Prepare holders for links in the sparse matrix
n_links = sum( target_indices ~= -1 );

% Put it in the adjacency matrix
index = 1;
for j = 1 : numel(target_indices)

% If we did not find a proper target to link, we skip
if target_indices(j) == -1
continue
end

% The source line number in the adjacency matrix
row_indices{i}(index) = current_slice_index + j;

% The target column number in the adjacency matrix
column_indices{i}(index) = current_slice_index + n_cells(i) + target_indices(j);

index = index + 1;

end

current_slice_index = current_slice_index + n_cells(i);

end

row_index = vertcat(row_indices{:});
column_index = vertcat(column_indices{:});
n_total_cells = sum(n_cells);

if debug
end

A = sparse(row_index, column_index, link_flag, n_total_cells, n_total_cells);

if debug
fprintf('Done.\n')
end

%% Gap closing

if debug
fprintf('Gap-closing:\n')
end

current_slice_index = 0;
for i = 1 : n_slices-2

% Try to find a target in the frames following, starting at i+2, and
% parsing over the target that are not part in a link already.

current_target_slice_index = current_slice_index + n_cells(i) + n_cells(i+1);

for j = i + 2 : min(i +  max_gap_closing, n_slices)

source = points{i}(unmatched_sources{i}, :);
target = points{j}(unmatched_targets{j}, :);

if isempty(source) || isempty(target)
current_target_slice_index = current_target_slice_index + n_cells(j);
continue
end

% Put it in the adjacency matrix
for k = 1 : numel(target_indices)

% If we did not find a proper target to link, we skip
if target_indices(k) == -1
continue
end

if debug
fprintf('Creating a link between point %d of frame %d and point %d of frame %d.\n', ...
unmatched_sources{i}(k), i, unmatched_targets{j}(target_indices(k)), j);
end

% The source line number in the adjacency matrix
row_index = current_slice_index + unmatched_sources{i}(k);
% The target column number in the adjacency matrix
column_index = current_target_slice_index + unmatched_targets{j}(target_indices(k));

A(row_index, column_index) = 1; %#ok<SPRIX>

end

new_links_target =  target_indices ~= -1 ;

current_target_slice_index = current_target_slice_index + n_cells(j);
end

current_slice_index = current_slice_index + n_cells(i);

end

if debug
fprintf('Done.\n')
end

%% Parse adjacency matrix to build tracks

if debug
fprintf('Building tracks:\n')
end

% Find columns full of 0s -> means this cell has no source
cells_without_source = [];
for i = 1 : size(A, 2)
if length(find(A(:,i))) == 0 %#ok<ISMT>
cells_without_source = [ cells_without_source ; i ]; %#ok<AGROW>
end
end

n_tracks = numel(cells_without_source);

AT = A';

for i = 1 : n_tracks

tmp_holder = NaN(n_total_cells, 1);

target = cells_without_source(i);
index = 1;
while ~isempty(target)
tmp_holder(index) = target;
target = find( AT(:, target), 1, 'first' );
index = index + 1;
end

adjacency_tracks{i} = tmp_holder ( ~isnan(tmp_holder) );
end

%% Reparse adjacency track index to have it right.
% The trouble with the previous track index is that the index in each
% track refers to the index in the adjacency matrix, not the point in
% the original array. We have to reparse it to put it right.

tracks = cell(n_tracks, 1);

for i = 1 : n_tracks

track = NaN(n_slices, 1);

for j = 1 : numel(adjacency_track)

% We must determine the frame this index belong to
tmp = cell_index;
frame_index = 1;
while tmp > 0
tmp = tmp - n_cells(frame_index);
frame_index = frame_index + 1;
end
frame_index = frame_index - 1;
in_frame_cell_index = tmp + n_cells(frame_index);

track(frame_index) = in_frame_cell_index;

end

tracks{i} = track;

end

end
```