%TrackMate Matlab graphing code
%Saumya Saurabh, Ph.D.
%saumya.saurabh@gmail.com
%Uses the cdfplot function

clear
% Specify the folder where the files live
%Ideally each folder will have data from a single strain/ mutant
myFolder = 'D:\Science\Projects\Explorations\CB15 flow\Images\SW Motility';
% Check to make sure that folder actually exists.  Warn user if it doesn't.
if ~isfolder(myFolder)
    errorMessage = sprintf('Error: The following folder does not exist:\n%s\nPlease specify a new folder.', myFolder);
    uiwait(warndlg(errorMessage));
    myFolder = uigetdir(); % Ask for a new one.
    if myFolder == 0
         % User clicked Cancel
         return;
    end
end


% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(myFolder, 'Copy_of_aln_pos2TrNew.xml'); % Change to whatever pattern you need.
analyzed_files = dir(filePattern);
mintracklength =19; % set this number to get the minimum num of traj points
newtracks = {};
for k = 1:length(analyzed_files)
    baseFileName = analyzed_files(k).name;
    file_path = fullfile(analyzed_files(k).folder, baseFileName);
    fprintf(1, 'Now reading %s\n', file_path);
    tracks = importTrackMateTracks(baseFileName);
    %remove tracks less than 20 frames
    newtracks{k,1} = tracks(cellfun('length',tracks)>mintracklength);     
end
tracks_consol = cat(1, newtracks{:});
%Use newtracks from here on to extract, analyze and plot
%Calculating the distances from origin
for m = 1:length(newtracks)
    for k = 1:length(newtracks{m,1})
        newtracks{m,1}{k,1}(:,5)= newtracks{m,1}{k,1}(:,2)-newtracks{m,1}{k,1}(1,2);
        newtracks{m,1}{k,1}(:,6)= newtracks{m,1}{k,1}(:,3)-newtracks{m,1}{k,1}(1,3);
    end
end

     
%Calculating the jumps in each track
for m = 1:length(newtracks)
    for k = 1:length(newtracks{m,1})
        newtracks{m,1}{k,1}(1,7) = 0;
        for p = 1:length(newtracks{m,1}{k,1})-1
            newtracks{m,1}{k,1}(p+1,7) = sqrt((newtracks{m,1}{k,1}(p+1,5)-newtracks{m,1}{k,1}(p,5)).^2 + (newtracks{m,1}{k,1}(p+1,6)-newtracks{m,1}{k,1}(p,6)).^2);
        end
    end
end
            
                
%collect all jumps from a sample in one matrix to use cdfplot
jump = [];
for m = 1:length(newtracks)
    for k = 1:length(newtracks{m,1})
        jump = [jump;newtracks{m,1}{k,1}(:,7)];
    end
end


jump=jump(jump~=0);
%subplot(2,1,2)


%Plotting the distances from origin
figure(221)
%subplot(2,1,1)
for m = 1:length(newtracks)
    for k = 1:length(newtracks{m,1})
        plot(newtracks{m,1}{k,1}(:,5),newtracks{m,1}{k,1}(:,6),'b-','LineWidth',2)
        hold on
    end
end
xlim([-50 50])
ylim([-50 50])
axis square



%MSD analyses using MSDAnalyzer
SPACE_UNITS = 'µm';
TIME_UNITS = 's';
ma = msdanalyzer(3, SPACE_UNITS, TIME_UNITS);
ma = ma.addAll(tracks_consol);
v = ma.getVelocities %#ok<NOPTS>
V = vertcat( v{:} );
V(:,4) = sqrt(V(:,2).^2+V(:,2).^2)/0.052;
h = histogram(V(:, 4),35,'Normalization','pdf');
h.FaceColor = [0 0.5 0.5];
h.EdgeColor = 'r';
set(gca,'Fontsize',30)
xlabel('Instantaneous Speed (\mum/s)') 
ylabel('#') 
title('\DeltaPAS')

ma = ma.computeVCorr;
ma.vcorr
ma.plotMeanVCorr


%Plotting the distances from origin
figure(221)
%subplot(2,1,1)
for m = 1:length(newtracks)
    for k = 1:length(newtracks{m,1})
        plot(newtracks{m,1}{k,1}(:,5),newtracks{m,1}{k,1}(:,6),'b-','LineWidth',2)
        hold on
    end
end
xlim([-50 50])
ylim([-50 50])
axis square


figure(225)
cdfplot(jump);

axis square

save('DPAS_analyzed.mat');



%plotting histograms of speed together
clear
close all
load('WT_analyzed.mat')
ma = ma.computeVCorr;
ma.vcorr
h = ma.plotMeanVCorr
h.Color = 'k';
hold on
clear
load('DPAS_analyzed.mat')
ma = ma.computeVCorr;
ma.vcorr
h = ma.plotMeanVCorr
h.Color = 'b';
set(gca,'Fontsize',30)
xlabel('Instantaneous Speed (\mum/s)') 
ylabel('Distance (\mum)') 



%plotting VACF together
clear
close all
load('WT_analyzed.mat')
h = histogram(V(:, 4),35);
h.FaceColor = [0 0.5 0.5];
h.EdgeColor = 'r';
hold on
clear
load('DPAS_analyzed.mat')
h = histogram(V(:, 4),30);
h.FaceColor = [0 0.9 0.9];
h.EdgeColor = 'b';
set(gca,'Fontsize',30)
xlabel('Instantaneous Speed (\mum/s)') 
ylabel('Distance (\mum)') 



%if multiple plots are needed consider running the above code block n times
%keep everything same just change XML directory and the variable jump to jump_WT, jump_DPAS etc
%cdfplot(jump_WT); hold on; cdfplot(jump_DPAS); hold on... etc. 

% %just for plotting swarms one by one
% clear all
% close all
% load('DPAS_analyzed.mat')
% for m = 1:length(newtracks)
%     for k = 1:length(newtracks{m,1})
%         plot(newtracks{m,1}{k,1}(:,5),newtracks{m,1}{k,1}(:,6),'b-','LineWidth',2)
%         hold on
%     end
% end
% xlim([-50 50])
% ylim([-50 50])
% set(gca,'Linewidth',3,'Fontsize', 30)
% xlabel('Distance (\mum)') 
% ylabel('Distance (\mum)') 
% title('\DeltaPAS')
% axis square
% 
% %just for plotting CDF together
% clear all
% close all
% load('DFliG_analyzed.mat')
% cdfplot(jump);
% hold on
% load('WT_analyzed.mat')
% cdfplot(jump);
% hold on
% load('DPAS_analyzed.mat')
% cdfplot(jump);
% hold on
% load('DCache_analyzed.mat')
% cdfplot(jump);


% %just for plotting velocity histograms together
clear all
close all
load('WT_analyzed.mat')
h = histogram(V(:, 4),32,'Normalization','probability');
%h = histogram(V(:, 4),35);
h.FaceColor = [0 0 1]; %Change color here
h.EdgeColor = 'b'; %Change color here
hold on
load('DPAS_analyzed.mat')
h = histogram(V(:, 4),32,'Normalization','probability');
h.FaceColor = [0.5 0.5 0.5]; %Change color here
h.EdgeColor = 'k'; %Change color here
hold on
set(gca,'Fontsize',30)
xlim([7,80])
ylim([0,0.035])
xlabel('Instantaneous Speed (\mum/s)') 
ylabel('#') 






