cell_postprocessing.m

function cell_postprocessing(props)

plot_dXdY_max = props.plot_dXdY_max;
plot_UV_max = props.plot_UV_max;
plot_cell_num_max = props.plot_cell_num_max;
d_acc_threshold = props.d_acc_threshold;
time_step = props.time_step;
case_string = props.case_string;

track_dir =  fullfile('cases',case_string,'tracks');
particle_dir =  fullfile('cases',case_string,'particles');
plot_dir = 'results';
mkdir(plot_dir)
load([track_dir,'/tracks.mat'])
num_frms = numel(dir(fullfile(particle_dir ,'*.mat')));
num_tracks = size(tracks_order,3);
close all
figure('position',[0 0 1000 1000],'color','w')

dX = [];
dY = [];
U = [];
V = [];
d_euclid = zeros(1,num_tracks);
d_accum = zeros(1,num_tracks);
xy_end = zeros(num_tracks,2);
c = 0;
flag = zeros(1,num_tracks);
motile_num = zeros(1,num_frms);
longest_length = 0;

for p = 1 : num_tracks
    
    % extract current Lagrangian track
    current_track = tracks_order(:,:,p);
    ind = current_track(:,2) ~= 0;
    longest_length = max(longest_length,sum(ind));
    current_track = current_track(ind,:);
    
    % subtract the original place
    X = current_track(:,2) - current_track(1,2);
    Y = current_track(:,3) - current_track(1,3);
    
    % instantaneous displacement over each frame
    U_p = X(2:end) - X(1:end-1);
    V_p = Y(2:end) - Y(1:end-1);
    d_accum_p = sum((U_p.^2 + V_p.^2).^0.5);
    
    if d_accum_p >= d_acc_threshold
        frames = find(ind);
        try
        for f = frames
            motile_num(f) = motile_num(f) + 1;
        end
        catch
            keyboard
        end
        c = c + 1;
        flag(p) = 1;
        d_accum(c) = d_accum_p;
        U = [U; U_p(:)];
        V = [V; V_p(:)];
    
        % integral displacement from the original place
        dX_p = X(end);
        dY_p = Y(end);
        xy_end(c,:) = [dX_p dY_p];
        d_euclid(c) = (dX_p^2 + dY_p^2)^0.5;
        dX = [dX; dX_p(:)];
        dY = [dY; dY_p(:)];
    end  
end

tracks_thresholded = tracks_order(:,:,find(flag));
save([track_dir,'/tracks_thresholded.mat'],'tracks_thresholded')

ind = 1:c;
d_euclid = d_euclid(ind);
d_accum = d_accum(ind); 
xy_end = xy_end(ind,:);
% Directness
D_p = d_euclid./d_accum; 
D_a = mean(D_p);

% Center of Mass
M_end = mean(xy_end,1);

% Forward Migration Index
y_FMI = mean(xy_end(:,2)./d_accum(:));
x_FMI = mean(xy_end(:,1)./d_accum(:));

time_axis = (0:num_frms-1)  * time_step/60;

cell_num = zeros(num_frms,1);
file_list = dir(fullfile(particle_dir,'*.mat'));
for t = 1: num_frms
    file_name = fullfile(particle_dir,file_list(t).name);
    load(file_name);    
    cell_num(t) = length(particles.pos);  
end

results = [];
results.time_step = time_step;
results.case_string = case_string;
results.d_acc_threshold = d_acc_threshold;

% results.mean_int_threshold = basic_t;
% results.min_peak_int_threshold = peak_min;
% results.min_area = area_min;

results.U = U;
results.V = V;
results.dX = dX;
results.dY = dY;
results.cell_num = cell_num;
results.motile_num = motile_num;

results.d_accum = d_accum; % added by Brian
results.D_a = D_a;
results.M_end = M_end;
results.FMI = [x_FMI,y_FMI];

save(fullfile(plot_dir,case_string),'results')
%==========================================================================
subplot(339)
plot(0,0)
text(0,8,['Case: ',strrep(case_string,'_',' ')])
text(0,6,sprintf('Averaged directness = %.2f',D_a))
text(0,4,sprintf('Averaged center of Mass = (%.2f, %.2f)',M_end(1),M_end(2)))
text(0,2,sprintf('Forward Migration Index = (%.2f, %.2f)',x_FMI,y_FMI))
axis equal off
xlim([0 10])
ylim([0 10])

subplot(338)
cell_num_s = smooth(cell_num,0.2);
motile_num_s = smooth(motile_num,0.2);
plot(time_axis,cell_num_s,'O','linewidth',2,'markersize',5); hold on;
plot(time_axis,motile_num_s,'O','linewidth',2,'markersize',5); hold on;
xlim([0 time_axis(end)])
ylim([0 plot_cell_num_max])
title('Cell counts over time')
grid on;
xlabel('t (hours)')
ylabel('Cell counts')
legend('Total cell number','Motile cell number','location','southoutside')

%==========================================================================


% histogram
% dXdY_max = round(max(abs([dX(:);dY(:)])),-1);
% bin_size = 10;
% bins_half = (0:bin_size:dXdY_max) + bin_size/2;
% bins = [bins_half - (dXdY_max + bin_size),bins_half];
% [xx,yy] = meshgrid(edges,edges);
edges = - plot_dXdY_max :10: plot_dXdY_max;
N = hist3([dX,dY],'Ctrs',{edges,edges})';
N_x = sum(N,1);
N_y = sum(N,2);

subplot(331)
bar(edges,N_x/sum(N_x),'EdgeColor','none')
xlim([-1 1]*plot_dXdY_max)
title('Histogram of (X_f-X_0)')
xlabel('X_f-X_0 (pixel)')
grid on;

subplot(333)
bar(edges,N_y/sum(N_y),'EdgeColor','none')
xlim([-1 1]*plot_dXdY_max)
title('Histogram of (Y_f-Y_0)')
xlabel('Y_f-Y_0 (pixel)')
grid on;

subplot(332)
% color by tracks
colors_1 = distinguishable_colors(c);

% color by time - specific colormap
colors_2 = flipud(jet(longest_length));

% color by time - manual colormap
cc = (longest_length:-1:1)/longest_length;
colors_3 = [cc(:),cc(:),ones(longest_length,1)];

colors = colors_2;

p_set = find(flag);
for p = 1:20 %c
    
    % extract current Lagrangian track
    current_track = tracks_order(:,:,p_set(p));
    ind = current_track(:,2) ~= 0;
    current_track = current_track(ind,:);
    
    % subtract the original place
    X = current_track(:,2) - current_track(1,2);
    Y = current_track(:,3) - current_track(1,3);    

    % color by time - specific colormap
    for t = 2:length(X)
        plot([X(t) X(t-1)],[Y(t) Y(t-1)],'color',[colors(t,:),0.5],'linewidth',2,'linesmoothing','on'); hold on;
    end    
    plot(X(end),Y(end),'Ob',...
        'markerfacecolor','b','linewidth',2,'linesmoothing','on','markersize',2); hold on;
    
% color by tracks
%     plot(X,Y,'color',colors(p,:),'linewidth',1,'linesmoothing','on'); hold on;
%     plot(X(end),Y(end),'O','color',colors(p,:),...
%         'markerfacecolor',colors(p,:),'linewidth',2,'linesmoothing','on','markersize',2); hold on;
end
plot([0 0],[-1 1]*plot_dXdY_max,'--r','linewidth',2);  hold on;
plot([-1 1]*plot_dXdY_max,[0 0],'--r','linewidth',2);
xlabel('X - X_0 (pixel)')
ylabel('Y - Y_0 (pixel)')
axis equal
xlim([-1 1]*plot_dXdY_max)
ylim([-1 1]*plot_dXdY_max)
title('Cell movement from original place')


% UV_max = round(max(abs([U(:);V(:)])),-1);
% bin_size = 1;
% bins_half = (0:bin_size:UV_max) + bin_size/2;
% bins = [bins_half - (UV_max + bin_size),bins_half];
% [xx,yy] = meshgrid(bins,bins);
edges = - plot_UV_max : plot_UV_max;
N = hist3([U,V],'Ctrs',{edges,edges})';
N_x = sum(N,1);
N_y = sum(N,2);

n = props.roseBins;
[theta,rho] = cart2pol(U,V);
r = (0:n)'/n * plot_UV_max * sqrt(2);
t = pi*(-n:n)/n;
X = r*cos(t);
Y = r*sin(t);
N = hist3([theta,rho],'edges',{t,r+max(rho)/n})';

subplot(335)
% h = pcolor(xx-bin_size/2,yy-bin_size/2,N); hold on
h = pcolor(X,Y,N); colorbar; axis equal tight; hold on
%set(h, 'EdgeColor', 'none')
%shading interp
plot([0 0],[-200 200],'--r','linewidth',2);  hold on;
plot([-200 200],[0 0],'--r','linewidth',2);
xlabel('U (pixel/frame)')
ylabel('V (pixel/frame)')
axis equal
xlim([-1 1]*plot_UV_max)
ylim([-1 1]*plot_UV_max)
title({'2D histogram',' of movement per frame'})
hold on
colorbar

subplot(334)
bar(edges,N_x/sum(N_x),'EdgeColor','none')
xlim([-1 1]*plot_UV_max)
title('Histogram of U')
xlabel('U (pixel/frame)')
grid on;

subplot(336)
bar(edges,N_y/sum(N_y),'EdgeColor','none')
xlim([-1 1]*plot_UV_max)
title('Histogram of V')
xlabel('V (pixel/frame)')
grid on;
print(fullfile(plot_dir,case_string),'-dtiff','-r300')
%==========================================================================