%This script generates the images for figure 7E 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=[7 15 31 47 63]
%% -----set up------
folder=['../data/flow_phantom/d/UIV/2cm/'];
name=['/phantom_d_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=[21:50 1:20];
        centre=[1.9,1.0];
    elseif delta==15
        indices=[2:50 1];
        centre=[1.9,1.0];
    elseif delta==31
        indices=[44:50 1:43];
        centre=[1.9,1.0];
    elseif delta==47
        indices=[18:50 1:17];
        centre=[1.9,1.0];
    elseif delta==63
        indices=[24:50 1:23];
        centre=[1.9,1.0];
    end

%--sort out the time steps--
lines=128;
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/d/flow_rate/phantom_d_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=415.51./mean(flowT);
flowTMax=flowTMax.*scaleFactor;
flowTMin=flowTMin.*scaleFactor;
flowT=flowT.*scaleFactor;

%---flow rate UIV---
[flow2, Flow,upperEdge,lowerEdge]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 380, 430, 680, 730);
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
Tindices=[length(flowT)-180:length(flowT) 1:length(flowT)-181];
flowT2 = interp1([timeT], [flowT(Tindices)],times(1:end-1));
meanFlowerror=mean(abs(flow(1:end-1)-flowT2'))
percentFlowerror=meanFlowerror./(max(flow)).*100
%% -----velocity profile---------
clear vprofile
yfrom=round(mean(mean(upperEdge./(iastep2)))); yto=round(mean(mean(lowerEdge./(iastep2))));
n=1
for t=indices
    vprofile(:,n)=mean(-VxC([1:22 36:63],yfrom:yto,t),1);
    n=n+1;
end
r=([1:length(vprofile(:,t))]).*iastep2.*ypix.*-1;
r=(r-r(end)/2)*1000;

%% ----womersley profile----
figure;
for t=1:50
    R(t)=ypix*mean(lowerEdge(:,t)-upperEdge(:,t))/2;
end
Rmean=mean(R);
flow_f=flowT(Tindices);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
vmean=flow_f(pos)'./(pi.*R(indices).^2);
vmean2=vmean./(60*1000*1000);
nu=1e-6;
vel= womer_Vmean(timeT(pos), 20, vmean2, Rmean,'j',nu);

flow_f(pos)'./(pi.*Rmean.^2);
vmean2=vmean./(60*1000*1000);
vel3= womer_Vmean(times, 20, vmean2, Rmean,nu);

vprofile_average=(vprofile(ceil(length(vprofile)/2):length(vprofile)-1,:)+vprofile(fliplr(1:floor(length(vprofile)/2)),:))/2;

figure;  plot([1:50], real(vel3(:,1)),'k');
hold on; plot([1:50], vprofile(round(length(vprofile)/2), :),'r');
legend('Womersley at R=0', 'UIV at R=0');
figfilename=['phantom_d_Delta' int2str(delta) '_WomersleyComparison_R=0b.png'];
print( '-painters', '-dpng', '-r600', figfilename);

womerpro=real(vel3);
UIVpro=interp1(r(fliplr(1:floor(length(vprofile)/2))),vprofile_average(:,:),[0:0.01:1].*R(indices(t))*1000,'spline','extrap');
meanWomerror=mean(mean(abs(womerpro-UIVpro')))
percenterror=meanWomerror./(max(max(womerpro))).*100
percenterrorAxis=abs(womerpro(:,1)-UIVpro(1,:)')./abs(womerpro(:,1));

%% -----1 velocity profile---------
yfrom=1; yto=241;
idx=21; 
for t=indices(idx)
    vmag=(VxC(:,:,t).^2+VyC(:,:,t).^2).^0.5;
    vmag=vmag.*sign(VxC(:,:,t));
    vprofile=mean(-vmag([11:45],yfrom:yto),1);
end
figure(102);
r=([1:length(vprofile)]).*iastep2.*ypix.*1;
Upper=round(mean(upperEdge(:,idx))/iastep2);
Lower=round(mean(lowerEdge(:,idx))/iastep2);
r=(r-ypix*iastep2*(Upper+Lower)/2)*1000;
r=r(Upper:Lower);

hold on;
if delta==7
    plot(r,vprofile(Upper:Lower),'Color',[0.55 0 0.55],'LineWidth',0.5);
elseif delta==15
    plot(r,vprofile(Upper:Lower),'m','LineWidth',0.5);
elseif delta==31
    plot(r,vprofile(Upper:Lower),'r','LineWidth',0.5);
elseif delta == 47
    plot(r,vprofile(Upper:Lower),'Color',[1 0.5 0],'LineWidth',0.5);
elseif delta==63
    plot(r,vprofile(Upper:Lower),'g','LineWidth',0.5);
elseif delta==79
    plot(r,vprofile(Upper:Lower),'c','LineWidth',0.5);
elseif delta==95
    plot(r,vprofile(Upper:Lower),'Color',[0.55 0 0.55],'LineWidth',0.5);
end


%% ------write out errors file
filename=['phantom_d_Delta' int2str(delta) 'errors.txt'];

fid=fopen(filename,'w');
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanFlowerror,percentFlowerror);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanWomerror,percenterror);
fclose('all');
end

%% ---conventional----
%% -----set up------
name=['../data/flow_phantom/d/UIV/2cm/conventional/phantom_d_128lines_conventional_20000depth7000000freq_out']
load([name '.mat']);

lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(lines, depth, Vx, Vy, 0,'water',period, split);
indices=[15:50 1:14];
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,upperEdge,lowerEdge]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 380, 430, 680, 730);
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 2600]);
xlabel('time /s','FontSize',8);
ylabel('flow rate /ml min^{-1}','FontSize',8);
strings={'\Delta t=7','\Delta t=15','\Delta t=31','\Delta t=47','\Delta t=63','normal','flow meter'};
legend(strings,'Position',[0.73 0.66 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_d_flow_rate_comparison.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%% -----velocity profile---------
clear vprofile
yfrom=round(mean(mean(upperEdge./(iastep2)))); yto=round(mean(mean(lowerEdge./(iastep2))));
n=1
for t=indices
    vprofile(:,n)=mean(-VxC([1:22 36:63],yfrom:yto,t),1);
    n=n+1;
end
r=([1:length(vprofile(:,t))]).*iastep2.*ypix.*-1;
r=(r-r(end)/2)*1000;

%% ----womersley profile----
figure;
for t=1:50
    R(t)=ypix*mean(lowerEdge(:,t)-upperEdge(:,t))/2;
end
Rmean=mean(R);
flow_f=flowT(Tindices);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
vmean=flow_f(pos)'./(pi.*R(indices).^2);
vmean2=vmean./(60*1000*1000);
vel= womer_Vmean(timeT(pos), 20, vmean2, Rmean,'j',nu);

flow_f(pos)'./(pi.*Rmean.^2);
vmean2=vmean./(60*1000*1000);
vel3= womer_Vmean(times, 20, vmean2, Rmean,'j',nu); %timeT(pos) vmean2(pos)

vprofile_average=(vprofile(ceil(length(vprofile)/2):length(vprofile)-1,:)+vprofile(fliplr(1:floor(length(vprofile)/2)),:))/2;

figure;  plot([1:50], real(vel3(:,1)),'k');
hold on; plot([1:50], vprofile(round(length(vprofile)/2), :),'r');
legend('Womersley at R=0', 'UIV at R=0');
figfilename=['phantom_d_conventional_WomersleyComparison_R=0b.png'];
print( '-painters', '-dpng', '-r600', figfilename);

womerpro=real(vel3);
UIVpro=interp1(r(fliplr(1:floor(length(vprofile)/2))),vprofile_average(:,:),[0:0.01:1].*R(indices(t))*1000,'spline','extrap');
meanWomerror=mean(mean(abs(womerpro-UIVpro')))
percenterror=meanWomerror./(max(max(womerpro))).*100
percenterrorAxis=abs(womerpro(:,1)-UIVpro(1,:)')./abs(womerpro(:,1));

%% -----1 velocity profile---------
yfrom=1; yto=241;
for t=indices(idx)
    vmag=(VxC(:,:,t).^2+VyC(:,:,t).^2).^0.5;
    vmag=vmag.*sign(VxC(:,:,t));
    vprofile=mean(-vmag([11:45],yfrom:yto),1);
end
figure(102);
r=([1:length(vprofile)]).*iastep2.*ypix.*1;
Upper=round(mean(upperEdge(:,idx))/iastep2);
Lower=round(mean(lowerEdge(:,idx))/iastep2);
r=(r-ypix*iastep2*(Upper+Lower)/2)*1000;
r=r(Upper:Lower);
xlim([-3 3]);
ylim([0 2]);
hold on;
plot(r,vprofile(Upper:Lower),'b','LineWidth',0.5);
set(gca,'Position', [0.14 0.2 0.53 0.7], 'FontSize',8);
title(['time = ' num2str(idx/50.*times(end),2) ' s'],'FontSize',8);
xlabel('radial position /mm','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
strings={'\Delta t=7','\Delta t=15','\Delta t=31','\Delta t=47','\Delta t=63','normal'};
legend(strings,'Position',[0.75 0.5 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_d_peak_vel_profile_comparison.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%% ------write out errors file
filename=['phantom_d_conventional_errors.txt'];
fid=fopen(filename,'w');

formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanFlowerror,percentFlowerror);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanWomerror,percenterror);
fclose('all');