%Code to plot vectors calculated by FIDVC (Bar-Kochba E., Toyjanova J., Andrews E., Kim K., Franck C. (2014) A fast 
% iterative digital volume correlation algorithm for large deformations. 
% Experimental Mechanics. doi: 10.1007/s11340-014-9874-2)
%
%Stadler et al. (2022) "A circular zone of attachment to the
%extracellular matrix provides directionality to the motility of Toxoplasma gondii in 3D"


%close all
load('PLF101003_1_60_inc.mat') %load file from FIDVC

%Pick timepoint comparsion and z-slice
t=7;
z_slice=5;

u1a = (u{t,1}{1,1});
u2a = (u{t,1}{1,2});
u3a = (u{t,1}{1,3});
u4a = (u{t,1}{1,4});


%scale u1, u2, u3 to be pixels
u1 = u1a*0.217;
u2 = u2a*0.217;
u3 = u3a*0.5;

%% scale u4 to be in pixels using u1, u2, and u3
u4_scale = zeros(49, 49, 7);

i = 1;
k = 1;
j = 1;
tt = t+1;
x=50;
y=50;
z=8;

while t < tt;
    while j < z;
        while k < y;
            while i < x;
                u4_scale(i,k,j) = sqrt(((u1(i,k,j))^2) + ((u2(i,k,j))^2) + ((u3(i,k,j))^2));
                u1_u2_mag(i,k,j) = sqrt(((u1(i,k,j))^2) + ((u2(i,k,j))^2));
                i = i + 1;
            end
            k = k + 1;
            i = 1;
        end
        j = j + 1;
        k = 1;
    end
    t = t + 1;
    j = 1;
    
end

%% Create arrays for above background 3D
ii = 1;
kk = 1;
jj = 1;


while jj < z;
    while kk < y;
        while ii < x;
            if u4_scale(ii,kk,jj) > 0.07; %calcuate based on no parasite moving during file
                u1_AboveBack(ii,kk,jj) = u1(ii,kk,jj);
                u2_AboveBack(ii,kk,jj) = u2(ii,kk,jj);
                u3_AboveBack(ii,kk,jj) = u3(ii,kk,jj);
                ii = ii + 1;
            else
                u1_AboveBack(ii,kk,jj) = 0;
                u2_AboveBack(ii,kk,jj) = 0;
                u3_AboveBack(ii,kk,jj) = 0;
                ii = ii + 1;
            end
        end
        kk = kk + 1;
        ii = 1;
    end
    jj = jj + 1;
    kk = 1;
end


%% Create arrays for above background 2D
ii2 = 1;
kk2 = 1;
jj2 = 1;


while jj2 < z;
    while kk2 < y;
        while ii2 < x;
            if u1_u2_mag(ii2,kk2,jj2) > 0.03; %calcuate based on no parasite moving during file
                u1_AboveBack2(ii2,kk2,jj2) = u1(ii2,kk2,jj2);
                u2_AboveBack2(ii2,kk2,jj2) = u2(ii2,kk2,jj2);
                u3_AboveBack2(ii2,kk2,jj2) = u3(ii2,kk2,jj2);
                ii2 = ii2 + 1;
            else
                u1_AboveBack2(ii2,kk2,jj2) = 0;
                u2_AboveBack2(ii2,kk2,jj2) = 0;
                u3_AboveBack2(ii2,kk2,jj2) = 0;
                ii2 = ii2 + 1;
            end
        end
        kk2 = kk2 + 1;
        ii2 = 1;
    end
    jj2 = jj2 + 1;
    kk2 = 1;
end

%Make arrays at one z slice for all arrows (use for 2D)
u1z = u1(:,:,z_slice);
u2z = u2(:,:,z_slice);
u3z = u3(:,:,z_slice);

%Make arrays at one z slice for above background arrows (use for 2D)
u1z_AB = u1_AboveBack2(:,:,z_slice);
u2z_AB = u2_AboveBack2(:,:,z_slice);
u3z_AB = u3_AboveBack2(:,:,z_slice);


%%PLOT
f1 = figure;
f2 = figure;
f3 = figure;
f4 = figure;


arrow_length = 1; %arrow_length is a scale factor
scale = 1; %scales the size of the arrows
arrow_thickness = 0.5; %how thick the arrows are
head_size = 1 %arrow head size

hold on
%Plot above background arrows only 3D
figure(f1)
[Xgrid, Ygrid, Zgrid]=meshgrid(1:1.736:85.064, 1:1.736:85.064, 1:3.5:23.5);
quiver3(Xgrid,Ygrid,Zgrid,u1_AboveBack, u2_AboveBack, u3_AboveBack, 'AutoScaleFactor', arrow_length, 'Linewidth', arrow_thickness,'MaxHeadSize',head_size)
axis equal

%Plot all arrows 3D
figure(f2)
[Xgrid, Ygrid, Zgrid]=meshgrid(1:1.736:85.064, 1:1.736:85.064, 1:3.5:23.5);
quiver3(Xgrid,Ygrid,Zgrid,u1, u2, u3,'AutoScaleFactor', arrow_length, 'Linewidth', arrow_thickness,'MaxHeadSize',head_size)
axis equal

%plot all 2D arrows at chosen z-slice
figure(f3)
[Xgrid, Ygrid]=meshgrid(0.868:1.736:84.196);
quiver(Xgrid, Ygrid, u1z, u2z, 'AutoScaleFactor', 1, 'Linewidth', 1.5,'MaxHeadSize',0.2)
axis square
axis([0,85.932,0,85.932]);

%Plot above background 2D arrows at chosen z-slice
figure(f4)
[Xgrid, Ygrid]=meshgrid(0.868:1.736:84.196);
quiver(Xgrid, Ygrid, u1z_AB,u2z_AB, 'AutoScaleFactor', 1, 'Linewidth', 1.5,'MaxHeadSize',0.2)
axis square
axis([0,85.932,0,85.932]);


