clear all, close all

% Path to folder containing "import_trackmate_xml.m"
addpath %yourpath

% Input minimum number of timepoints required for tracks to be included in analysis
min_timepoints = 15;

% Input minima and maxima for total distance traveled and directionality of
% single-cell tracks. Only tracks that meet these criteria will be included
% in filtered results.
min_travel = 0;
max_travel = 500;
min_angle = 0;
max_angle = 360;

myDir = pwd;
myFiles = dir(fullfile(myDir,'*.xml'));


tables = cell(4,length(myFiles));

w = waitbar(0, 'One moment...');

for l = 1:length(myFiles)
  baseFileName = myFiles(l).name;
  fullFileName = fullfile(myDir, baseFileName);
  
  waitbar(l/length(myFiles), w, ['Now reading ' baseFileName]);
  
  % Use function to read xmltracks into MATLAB
  MTD = import_trackmate_xml(baseFileName);
  for i =1:height(MTD)
     
    % Set up temporary tables for each variable
    Timetemp1 = MTD.Time(i,1);
    Xtemp1 = MTD.X(i,1);
    Ytemp1 = MTD.Y(i,1); 
    Timetemp = Timetemp1{1,1};
    Xtemp = Xtemp1{1,1};
    Ytemp = Ytemp1{1,1};
    
    % Remove all datapoints where X & Y positions are both whole
    % numbers (TrackMate artifacts when it loses a cell trace)
    Lremainders = rem(Xtemp,1);       
    Mremainders = rem(Ytemp,1);
    Particleweird = Lremainders == 0 & Mremainders == 0;
    Timetemp(:,Particleweird) = [];
    Xtemp(:,Particleweird) = [];
    Ytemp(:,Particleweird) = [];
    
    % Skip all tracks with too few timepoints tracked

    if numel(Timetemp) >= min_timepoints

            
            
            % Calculate changes in position
            DeltaPostemp = nan(numel(Xtemp),1);
            for k = 2:numel(Xtemp)
                DeltaPostemp(k) = sqrt((Xtemp(k) - Xtemp(k-1))^2 + (Ytemp(k) - Ytemp(k-1))^2)/(Timetemp(k) - Timetemp(k-1));
                
                dX=Xtemp(k) - Xtemp(k-1);
                dY=Ytemp(k) - Ytemp(k-1);
                if dX>0&&dY>0
                    Angletemp(k-1)=atan(abs(dY/dX))*180/pi;
                elseif dX<0&&dY>0
                    Angletemp(k-1)=atan(abs(dX/dY))*180/pi+90;
                elseif dX<0&&dY<0
                    Angletemp(k-1)=atan(abs(dY/dX))*180/pi+180;
                elseif dX>0&&dY<0
                    Angletemp(k-1) = atan(abs(dX/dY))*180/pi+270;
                end
            end
            
            % Calculate angle of travel between all sets of two consecutive
            % datapoints
            DX = Xtemp(numel(Timetemp)) - Xtemp(1);
            DY = Ytemp(numel(Timetemp)) - Ytemp(1);
            if dX>0&&dY>0
                Angletemp_total = atan(abs(dY/dX))*180/pi;
            elseif dX<0&&dY>0
                Angletemp_total = atan(abs(dX/dY))*180/pi+90;
            elseif dX<0&&dY<0
                Angletemp_total = atan(abs(dY/dX))*180/pi+180;
            elseif dX>0&&dY<0
                Angletemp_total = atan(abs(dX/dY))*180/pi+270;
            end
            
            % Output cleaned up tracking data and summaries for each cell
            Tracks(i).ParticleID = MTD.ParticleID(i);
            Tracks(i).Time = Timetemp;
            Tracks(i).X = Xtemp;
            Tracks(i).Y = Ytemp;
            Tracks(i).MSD = DeltaPostemp';
            Tracks(i).SingleCellMeanMotility = mean(DeltaPostemp, 'omitnan');
            Tracks(i).SingleCellSTDMotility = std(DeltaPostemp, 'omitnan');
            Tracks(i).SingleCellTotalDistance = sum(DeltaPostemp, 'omitnan');
            Tracks(i).Angle = Angletemp;
            Tracks(i).SingleCellMeanAngle = mean(Angletemp);
            Tracks(i).SingleCellTotalAngle = Angletemp_total;
            
            if min_angle < Angletemp_total && Angletemp_total < max_angle && min_travel < Tracks(i).SingleCellTotalDistance && Tracks(i).SingleCellTotalDistance < max_travel
                filterTracks(i).ParticleID = MTD.ParticleID(i);
                filterTracks(i).Time = Timetemp;
                filterTracks(i).X = Xtemp;
                filterTracks(i).Y = Ytemp;
                filterTracks(i).MSD = DeltaPostemp';
                filterTracks(i).SingleCellMeanMotility = mean(DeltaPostemp, 'omitnan');
                filterTracks(i).SingleCellSTDMotility = std(DeltaPostemp, 'omitnan');
                filterTracks(i).SingleCellTotalDistance = sum(DeltaPostemp, 'omitnan');
                filterTracks(i).Angle = Angletemp;
                filterTracks(i).SingleCellMeanAngle = mean(Angletemp);
                filterTracks(i).SingleCellTotalAngle = Angletemp_total;
                
            else
                filterTracks(i).ParticleID = MTD.ParticleID(i);
                filterTracks(i).Time = Timetemp;
                filterTracks(i).X = Xtemp;
                filterTracks(i).Y = Ytemp;
                filterTracks(i).MSD = DeltaPostemp';
                filterTracks(i).SingleCellMeanMotility = [];
                filterTracks(i).SingleCellTotalDistance = [];
                filterTracks(i).Angle = [];
                filterTracks(i).SingleCellMeanAngle = [];
                filterTracks(i).SingleCellTotalAngle = [];
                
            end
            
    else
	    Tracks(i).ParticleID = MTD.ParticleID(i);
        Tracks(i).Time = Timetemp;
        Tracks(i).X = Xtemp;
        Tracks(i).Y = Ytemp;
        Tracks(i).MSD = DeltaPostemp';
        Tracks(i).SingleCellMeanMotility = [];
	    Tracks(i).SingleCellTotalDistance = [];
        Tracks(i).Angle = [];
        Tracks(i).SingleCellMeanAngle = [];
        Tracks(i).SingleCellTotalAngle = [];
        filterTracks(i).ParticleID = MTD.ParticleID(i);
        filterTracks(i).Time = Timetemp;
        filterTracks(i).X = Xtemp;
        filterTracks(i).Y = Ytemp;
        filterTracks(i).MSD = DeltaPostemp';
        filterTracks(i).SingleCellMeanMotility = [];
        filterTracks(i).SingleCellTotalDistance = [];
        filterTracks(i).Angle = [];
        filterTracks(i).SingleCellMeanAngle = [];
        filterTracks(i).SingleCellTotalAngle = [];

    end
    
     % Save data for imaging site in output cell
tables{1,l} = baseFileName(1:5);
tables{2,l} = Tracks;
tables{3,l} = mean([Tracks.SingleCellMeanMotility], 'omitnan');
tables{4,l} = std([Tracks.SingleCellMeanMotility], 'omitnan');
tables{5,l} = mean([Tracks.SingleCellTotalDistance], 'omitnan');
tables{6,l} = std([Tracks.SingleCellTotalDistance], 'omitnan');
tables{7,l} = mean([Tracks.SingleCellTotalAngle], 'omitnan');
tables{8,l} = std([Tracks.SingleCellTotalAngle], 'omitnan');
tables{9,l} = mean([filterTracks.SingleCellMeanMotility], 'omitnan');
tables{10,l} = std([filterTracks.SingleCellMeanMotility], 'omitnan');
tables{11,l} = mean([filterTracks.SingleCellTotalDistance], 'omitnan');
tables{12,l} = std([filterTracks.SingleCellTotalDistance], 'omitnan');
tables{13,l} = mean([filterTracks.SingleCellTotalAngle], 'omitnan');
tables{14,l} = std([filterTracks.SingleCellTotalAngle], 'omitnan');


 end
 
clear Tracks 
clear filterTracks
  
end

close(w)

label = {'File name'; 'All tracks'; 'Bulk mean motility'; 'Bulk stdev motility'; 'Bulk mean distance'; 'Bulk stdev distance';...
    'Bulk mean trajectory' ; 'Bulk stdev trajectory'; 'Filtered mean motility'; 'Filtered stdev motility'; 'Filtered mean distance'; 'Filtered stdev distance';...
    'Filtered mean trajectory' ; 'Filtered stdev trajectory';};
Tables = array2table(tables, 'RowNames', label);

save(['Tables',int2str(clock),'.mat'], 'Tables', '-v7.3')
