%This script generates the images for figure 8 and performs associated
%analysis.

%For full details see "Ultrasound Imaging Velocimetry with interleaved
%images for improved pulsatile arterial flow measurements: a new correction
%method, experimental and in vivo validation", J Royal Soc Interface, 2017

%K H Fraser, C Poelma 2016

for delta=[15 63]%for comparisson with deep
%% -----set up------
folder=['../data/flow_phantom/e/UIV/2cm/'];
name=['/phantom_e_128lines_' int2str(delta) 'delta_20000depth7000000freq_out']
load([folder int2str(delta) name '.mat']);
depth=0.02;
ypix=3.7425e-5/2;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, 1,'water',period, split);
    if delta==7
        indices=[19:50 1:18];
        centre=[1.9,1.3];
    elseif delta==15
        indices=[15:50 1:14];
        centre=[1.9,1.3];
    elseif delta==31
        indices=[27:50 1:26];
        centre=[1.9,1.3];
    elseif delta==47
        indices=[12:50 1:11];
        centre=[1.9,1.3];
    elseif delta==63
        indices=[11:50 1:10];
        centre=[1.9,1.3];
    end

%--sort out the time steps--
timestep=(lines*2+delta)*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;


%% ----flow rate-comparison----
% ----flow rate-transonics----
M=dlmread('../data/flow_phantom/e/flow_rate/phantom_e_flow.csv',',',[4 0 934 14]);
timeT=M(:,1);
timeT=timeT./1000;
flowT=M(:,12);
flowTMax=M(:,14);
flowTMin=M(:,15);

%scale by bucket flow rate to eliminate speed of sound error with
%Transonics
scaleFactor=634.2495./mean(flowT);
flowTMax=flowTMax.*scaleFactor;
flowTMin=flowTMin.*scaleFactor;
flowT=flowT.*scaleFactor;

%---flow rate UIV---
[flow2, Flow,upper,lower]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2,465, 530, 795, 855);
flow=flow2(indices);
figure(100); hold on;
if delta==7
    plot(times,flow,'Color',[0.55 0 0.55],'LineWidth',0.1);
elseif delta==15
    plot(times,flow,'m','LineWidth',0.1);
elseif delta==31
    plot(times,flow,'r','LineWidth',0.1);
elseif delta == 47
    plot(times,flow,'Color',[1 0.5 0],'LineWidth',0.1);
elseif delta==63
    plot(times,flow,'g','LineWidth',0.1);
elseif delta==79
    plot(times,flow,'c','LineWidth',0.1);
elseif delta==95
    plot(times,flow,'Color',[0.55 0 0.55],'LineWidth',0.1);
end

end

%% ---conventional----
%% -----set up------
name=['../data/flow_phantom/e/UIV/2cm/conventional/phantom_e_128lines_conventional_20000depth7000000freq_out']
load([name '.mat']);
depth=0.02;
ypix=3.7425e-5/2;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(lines, depth, Vx, Vy, 0,'water',period, split);
indices=[46:50 1:45];
centre=[1.9,1.3];

%--sort out the time steps--
lines=128;
timestep=(lines)*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;

%% ----flow rate-comparison----
%---flow rate UIV---
[flow2, Flow,upper,lower]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 465, 530, 795, 855);
flow=flow2(indices);
figure(100); hold on;
plot(times,flow,'b','LineWidth',0.1);

for delta=[7 31 63]
%% -----set up------
folder=['../data/flow_phantom/e/UIV/6cm/'];
name=['/phantom_e_128lines_' int2str(delta) 'delta_60000depth7000000freq_out']
load([folder int2str(delta) name '.mat']);
depth=0.06;
ypix=3.7425e-5/2;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, 1,'water',period, split);
    if delta==7
        indices=[27:50 1:26];
        centre=[1.9,1.3];
    elseif delta==15
        indices=[11:50 1:10];
        centre=[1.9,1.3];
    elseif delta==31
        indices=[40:50 1:39];
        centre=[1.9,1.3];
    elseif delta==47
        indices=[7:50 1:6];
        centre=[1.9,1.3];
    elseif delta==63
        indices=[15:50 1:14];
        centre=[1.9,1.3];
    end

%--sort out the time steps--
timestep=(lines*2+delta)*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;

%% ----flow rate-comparison----

%---flow rate UIV---
[flow2, Flow,upper,lower]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 465, 530, 795, 855);
flow=flow2(indices);
figure(100); hold on;
if delta==7
    plot(times,flow,'Color',[0.55 0 0.55],'linestyle','--','LineWidth',0.1);
elseif delta==15
    plot(times,flow,'m--','LineWidth',0.1);
elseif delta==31
    plot(times,flow,'r--','LineWidth',0.1);
elseif delta == 47
    plot(times,flow,'Color',[1 0.5 0],'linestyle','--','LineWidth',0.1);
elseif delta==63
    plot(times,flow,'g--','LineWidth',0.1);
elseif delta==79
    plot(times,flow,'c--','LineWidth',0.1);
elseif delta==95
    plot(times,flow,'Color',[0.55 0 0.55],'linestyle','--','LineWidth',0.1);
end

end


%% ---conventional----

%% -----set up------
name=['../data/flow_phantom/e/UIV/6cm/conventional/phantom_e_128lines_conventional_60000depth7000000freq_out']
load([name '.mat']);
depth=0.06;
ypix=3.7425e-5/2;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(lines, depth, Vx, Vy, 0,'water',period, split);
indices=[27:50 1:26];
centre=[1.9,1.3];

%--sort out the time steps--
timestep=(lines)*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;

%% ----flow rate-comparison----
%---flow rate UIV---
[flow2, Flow,upper,lower]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 290, 340, 570, 630);
flow=flow2(indices);
figure(100); hold on;
plot(times,flow,'b--','LineWidth',0.1);

set(gca,'Position',[0.19 0.2 0.6 0.76]);
Tindices=[length(flowT)-180:length(flowT) 1:length(flowT)-181];
hold on; plot(timeT, flowT(Tindices),'k','LineWidth',0.1);
hold on; plot(timeT,flowTMax(Tindices),'k:','LineWidth',0.1);
hold on; plot(timeT,flowTMin(Tindices),'k:','LineWidth',0.1);
%ylim([-100 1600]);
xlabel('time /s','FontSize',8);
ylabel('flow rate /ml min^{-1}','FontSize',8);
strings={'2 cm \Delta t=15 \tau','2 cm \Delta t=63 \tau','2 cm conventional','6 cm \Delta t=7 \tau','6 cm \Delta t=31 \tau','6 cm \Delta t=63 \tau','6 cm conventional','flow meter'};
legend(strings,'Position',[0.68 0.7 0.2 0.2],'FontSize',6);
legend('boxoff')
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [8.4 5]);
set(gcf, 'PaperPosition', [0 0 8.4 5]);
figfilename=['phantom_e_deep_flow_rate_comparison.png'];
print( '-painters', '-dpng', '-r600', figfilename);