%% process_tracks.m
% Written for Berro Lab at Yale University
% Last edited by Mike Lacy, January 2019

% This file contains two functions, process_tracks() and load_tracks().
% process_tracks() analyzes tracking data and groups datasets with
% processed features, it calls load_tracks() for each file.

% Some of the features are hard-coded based on the output format of PYME,
% but these are rather straightforward and could be adapted for other
% microscopy analysis files (we recorded movies in single image plane, so 
% motion metrics are calculated for 2D only).

% Example run command: >> sets = process_tracks({'fim1', 'act1'}, 100)
% It takes some time to run, several minutes or up to an hour for very
% large datasets. Run the script once and save the output sets as a .m file 
% for later use.

function sets = process_tracks(inputsets, frametime)
%% Inputs Summary
% inputsets{} is a cell array containing the names of the file sets to be 
% analyzed. Each file set has to be in a folder of the name given in 
% inputsets containing txt files generated from PYME tracking 
% (filenames must end in '_tracks.txt'), in the same parent directory as this script.
%
% *Be sure that all the tracks files you will analyze have used the same
% parameters in PYME (frametime, maxblinkframes, clumpradius, etc)
%
% frametime is the camera exposure time (in msec).

%% Outputs Summary
%sets() is a structure that contain all the tracks data, maintaining the
%organization by individual image files
% Its fields are given by:
% dirname: directory name of the file set.
% datasets: structure containing all the datasets pertaining to each sample.
% The fields of datasets() are:
%   filename: filename of the dataset.
%   tracks: structure containing information about the tracks.

% The following fields of tracks() come from the input files:
%   x,y,t: coordinates of each spot in each track, in nm and sec.
%   A: amplitude of each spot in each track, i.e. height of signal over background
%   sig: width of the gaussian fit to each spot in the track. FWHM ~ 2.355*sig
%   clumpIndex: clump Index of each track (i.e. track ID)
%   error_x, error_y: precision of localization in x and y, from PYME
%   plus any other spot features you choose to export from PYME

% Parameters calculated for each track include:
%   v: velocities of each spot in each track, for each step (nm/sec).
%   etedist: end-to-end-distance of each track (nm).
%   totdist: total displacement of each track (nm).
%   lifetime: lifetime of each track (sec).
%   offtime: array of all the offtime gaps recorded in each track (sec).
%   and others, feel free to inspect the code below or add new metrics as you like...

% NOTES 
% filenames/paths: Currently set up to run on Windows, where file paths
% use '\' instead of '/'.
% If you want to run on Mac, find all instances of '\' and replace with '/'
% - time data gets some weird float-point precision issues, 3/10 is =/= 0.3 

%% pre-settings
minlifetime = 0.3;% Minimum length of a track (in sec) allowed to include 
    % a track in the analysis. we already filtered for tracks of at least
    % points in PYME tracking, but we want a minimum time as well (ok if 
    % three points span 0.6 sec, but don't want a 0.3 sec event either).
maxallowedlifetime=40; % in sec.: to reject lifetimes longer than biologically reasonable (e.g. spot that persists the whole movie is probably junk or artifact). 
    % set to 0 if not using any limit.

dbstop if error;
tic;
%% Directories and filenames
% Start with the directory where the script is saved.
fullpath = mfilename('fullpath');
script_dir = strfind(fullpath,'\');
script_dir = fullpath(1:script_dir(end));

% Get all the filenames in each set's directory, 
% ignoring files who are not "tracks.txt"
sets=struct('dirname',inputsets);
N_sets=numel(sets);
for i=1:N_sets
    % Find all tracks.txt files in each dye directory.
    files=dir([script_dir,sets(i).dirname,'\*tracks.txt']);

    % Assign the filename to each datasets structure.
    allFilesNames={files.name};
    sets(i).datasets(1:numel(allFilesNames))=struct();
    [sets(i).datasets.filename]=deal(allFilesNames{:});
end

%% Run load_tracks() for each file and store results in the 'data' struct
% This loop is where the program spends *most* of its time: depends on how many tracks are in the datasets
for i=1:N_sets
    for j=1:numel(sets(i).datasets) 
        filename = strjoin({script_dir,sets(i).dirname,'\',sets(i).datasets(j).filename},'');
        data = load_tracks(filename, minlifetime, frametime);   
        % returns a struct with results of calculating track metrics
        % second input minlifetime is minimum allowed lifetime in seconds.
        
        if ~isempty(data)
            sets(i).datasets(j).tracks = data;
            % Option to reject tracks that are too long to be relevant 
            % (if known to be due to tracking errors or image artifacts)
            if(maxallowedlifetime)
                toolong = [sets(i).datasets(j).tracks.lifetime] > maxallowedlifetime;
                sets(i).datasets(j).tracks(toolong) = [];
            end
        end
    end
end
toc;
end % end of function process_tracks()

function data2 = load_tracks(filename, minlt, frametime)
% don't need to use this function, it is called by process_tracks().
% Import the tracks data from PYME tracking txt files.
fileID = fopen(filename);
vars_name=fgets(fileID);
vars_name=strrep(strsplit(vars_name),'#','');
vars_name(cellfun(@isempty,vars_name))=[];
N_vars=numel(vars_name);
tracks = textscan(fileID,repmat('%f ',1,N_vars));
fclose(fileID);

% Store the data in a structure.
data=struct();
N_datapoints=size(tracks{1},1);
for i=1:N_vars
    dat=num2cell(tracks{:,i});
    [data(1:N_datapoints).(vars_name{i})]=deal(dat{:});
end

last_track=max([data.clumpIndex]);

%% Group the data by their clump index. (they are pre-sorted by t)
data2=struct();
j=1;
for i=1:last_track
    ind=find([data.clumpIndex]==i);
    if max(ind)>0
        for k=1:N_vars
            data2(j).(vars_name{k})=[data(ind).(vars_name{k})];
        end
        data2(j).clumpIndex=i;
        
        % In some tracks, there are duplicate identifications of the same spot
        % In such cases, we choose the first entry in the .txt file and
        % discard the subsequent ones that appear in the same time frame.
        [~,ind,~]=unique([data2(j).t]);
        data2(j).x=data2(j).x(ind);
        data2(j).y=data2(j).y(ind);
        data2(j).t=data2(j).t(ind);
        
        % convert all time values from frames to seconds.
        data2(j).frametime = frametime;  %frame exposure time, in msec
        data2(j).t = data2(j).t .* (data2(j).frametime/1000);               
        % remember after this point, t will not be able to be compared by ==
        % because of floating point precision, i.e. 10 frames * 100msec is ~= 1.0 sec
        j=j+1;
    end
end

%% Metrics to calculate for each track

if ~isempty(fieldnames(data2))    
    % In cases with stringent tracking thresholds or sparse images 
    % sometimes there are no valid tracks from this file. Then data2 
    % would be an empty struct and would cause errors.
    
    initialt = min([data2(:).t]);
    for i=1:size(data2,2)
        data2(i).dx = diff(data2(i).x); %x, y, dx, dy are in units nm.
        data2(i).dy = diff(data2(i).y);
        data2(i).dt = diff(data2(i).t); % t and dt are in units of sec
        data2(i).trel = data2(i).t - initialt;  %adjusted relative time from first detected spot (start of illumination)
        data2(i).steps = sqrt(data2(i).dx.^2 + data2(i).dy.^2); %distance of successive steps, in nm
        data2(i).v = data2(i).steps ./ data2(i).dt; %instantaneous speed at each point, nm/sec
        data2(i).lifetime = data2(i).t(end)-data2(i).t(1);
        data2(i).start = data2(i).t(1);
        
        data2(i).meanint = mean(data2(i).A);
        
        data2(i).comx = mean(data2(i).x);   %center of mass, i.e. particle's time-averaged position.
        data2(i).comy = mean(data2(i).y);
        data2(i).xrel = data2(i).x - data2(i).comx; % position relative to center of mass
        data2(i).yrel = data2(i).y - data2(i).comy;
        
        data2(i).xete = data2(i).x(end) - data2(i).x(1);
        data2(i).yete = data2(i).y(end) - data2(i).y(1);
            %^the end-to-end vector is [xete yete]. *note the direction 
            % will not be informative if not aligned from cell edge vector*
        data2(i).etedist = sqrt((data2(i).x(end)-data2(i).x(1))^2 + (data2(i).y(end)-data2(i).y(1))^2);
        data2(i).totdisp = sum(sqrt((data2(i).dx).^2 + (data2(i).dy).^2));
        
        data2(i).offtime = data2(i).dt(abs(data2(i).dt - data2(i).frametime/1000) > 1e-6); 
            % any point where time between points is more than 1 frame.
            % to deal with floating point precision, allow relative error of 1e-6
    end
end
if (~isempty(fieldnames(data2)))
    data2([data2.lifetime]<minlt)=[]; % strip out any tracks shorter than the minimum lifetime
else
    data2 = []; 
    %if there are no tracks contained in data2, return an empty array instead of struct with empty values
end
end % end of function load_tracks()

