%This script generates the images for figure 10 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

delta=[31]
%% -----set up------
folder=['../data/flow_phantom/b/UIV/2cm/'];
name=['/phantom_b_128lines_' int2str(delta) 'delta_20000depth7000000freq_out']
load([folder int2str(delta) name '.mat']);
depth=0.02;
ypix=2.3653e-05;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_justscale(delta, depth, Vx, Vy, 1,'glyce',period, split);
if delta==31
    indices=[50 1:49];
    centre=[2.3,0.9];
    heightBox=0.1;
    widthBox=0.2;
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;
%---filter t---
%use for high pulse rate
VxC2(:,:,1)=VxC(:,:,50);
VxC2(:,:,2:51)=VxC(:,:,1:50);
VxC2(:,:,52)=VxC(:,:,1);
VyC2(:,:,1)=VyC(:,:,50);
VyC2(:,:,2:51)=VyC(:,:,1:50);
VyC2(:,:,52)=VyC(:,:,1);
[a,b,c]=size(VxC2);
for l=1:a
    for r=1:b
        VxC3(l,r,:)=medfilt1(squeeze(VxC2(l,r,:)),[3]);
        VyC3(l,r,:)=medfilt1(squeeze(VyC2(l,r,:)),[3]);
    end
end
VxC=VxC3(:,:,2:51);
VyC=VyC3(:,:,2:51);

%% ----flow rate-comparison----
%---flow rate UIV---
reg1=((centre(1)-widthBox/2)./3.8)*128;
reg2=((centre(1)+widthBox/2)./3.8)*128;
reg_lines=[reg1:reg2];
[flow2, Flow,upper,lower]=flowRatePulseA_f_reg(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 290, 340, 570, 630,reg_lines);
flow=flow2(indices);
figure(100); hold on;
plot(times,flow,'r--','LineWidth',0.1);

%% --simulate the Doppler spectrum--
%extract a box of data
figure;
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=imagesc(width,depth,log(1+abs(hilbert(squeeze(Iavg(iastep1:iastep1:end,iastep2:iastep2:end,1)))))');
axis equal;  axis tight;

[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-50:150;
intensity=zeros(length(velocities),50);
%simulate the spectrum within that box
T=1;
for t=[1:50]
    vmag=(VxC(:,:,t)'.^2+VyC(:,:,t)'.^2).^0.5;
    vmag=-1.*vmag.*sign(VxC(:,:,t)');
    pos=find([BW(1:end-1,1:end-1)]);
    vmag_small=vmag(pos);
    pos=find(isfinite(vmag_small));
    vmag_smaller=vmag_small(pos);
    v_nocorr(T)=mean(vmag_smaller)+std(vmag_smaller);
    T=T+1;
end

%spectra
figure(200);
hold on; plot(times, v_nocorr(indices),'r--');
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;

%% -----no acceleration correction
depth=0.02;
[VxC, VyC]=pulsePIVScaleAndCorrect_noacceleration(delta, depth, Vx, Vy, 1,'glyce',period, split);

%---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;
%---filter t---
VxC2(:,:,1)=VxC(:,:,50);
VxC2(:,:,2:51)=VxC(:,:,1:50);
VxC2(:,:,52)=VxC(:,:,1);
VyC2(:,:,1)=VyC(:,:,50);
VyC2(:,:,2:51)=VyC(:,:,1:50);
VyC2(:,:,52)=VyC(:,:,1);
[a,b,c]=size(VxC2);
for l=1:a
    for r=1:b
        VxC3(l,r,:)=medfilt1(squeeze(VxC2(l,r,:)),[3]);
        VyC3(l,r,:)=medfilt1(squeeze(VyC2(l,r,:)),[3]);
    end
end
VxC=VxC3(:,:,2:51);
VyC=VyC3(:,:,2:51);

%% ----flow rate-comparison----
%---flow rate UIV---
[flow2, Flow,upper,lower]=flowRatePulseA_f_reg(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 290, 340, 570, 630,reg_lines);
flow=flow2(indices);
figure(100); hold on;
plot(times,flow,'r:','LineWidth',0.1);
%% --simulate the Doppler spectrum--
%extract a box of data
figure;
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=imagesc(width,depth,log(1+abs(hilbert(squeeze(Iavg(iastep1:iastep1:end,iastep2:iastep2:end,1)))))');
axis equal;  axis tight;

[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-50:150;
intensity=zeros(length(velocities),50);
%simulate the spectrum within that box
T=1;
for t=[1:50]
    vmag=(VxC(:,:,t)'.^2+VyC(:,:,t)'.^2).^0.5;
    vmag=-1.*vmag.*sign(VxC(:,:,t)');
    pos=find([BW(1:end-1,1:end-1)]);
    vmag_small=vmag(pos);
    pos=find(isfinite(vmag_small));
    vmag_smaller=vmag_small(pos);
    v_vcorr(T)=mean(vmag_smaller)+std(vmag_smaller);
    T=T+1;
end

%spectra
figure(200);
hold on; plot(times, v_vcorr(indices),'r:');
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;

 %% ----complete correction
 depth=0.02;
 [VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, 1,'glyce',period, split);

%---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;
%---filter t---
VxC2(:,:,1)=VxC(:,:,50);
VxC2(:,:,2:51)=VxC(:,:,1:50);
VxC2(:,:,52)=VxC(:,:,1);
VyC2(:,:,1)=VyC(:,:,50);
VyC2(:,:,2:51)=VyC(:,:,1:50);
VyC2(:,:,52)=VyC(:,:,1);
[a,b,c]=size(VxC2);
for l=1:a
    for r=1:b
        VxC3(l,r,:)=medfilt1(squeeze(VxC2(l,r,:)),[3]);
        VyC3(l,r,:)=medfilt1(squeeze(VyC2(l,r,:)),[3]);
    end
end
VxC=VxC3(:,:,2:51);
VyC=VyC3(:,:,2:51);

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

M=dlmread('../data/flow_phantom/b/flow_rate/phantom_b_flow.csv',',',[4 0 288 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=365.48./mean(flowT);
flowTMax=flowTMax.*scaleFactor;
flowTMin=flowTMin.*scaleFactor;
flowT=flowT.*scaleFactor;

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

set(gca,'Position',[0.2 0.18 0.75 0.78]);
Tindices=[235:length(flowT) 1:234];
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([-200 1350]);
strings={'no corrections','with steady velocity correction','with total correction','flow meter'};
legend(strings,'Position',[0.6 0.78 0.2 0.14],'FontSize',6);
legend('boxoff')
xlabel('time /s','FontSize',8);
ylabel('flow rate /ml min^{-1}','FontSize',8);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [8.4 5]);
set(gcf, 'PaperPosition', [0 0 8.4 5]);
figfilename=['phantom_b_Delta' int2str(delta) '_flow_rate_acceleration_2.3.png'];
print( '-painters', '-dpng', '-r600', figfilename);
%% --simulate the Doppler spectrum--
%extract a box of data
figure;
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=imagesc(width,depth,log(1+abs(hilbert(squeeze(Iavg(iastep1:iastep1:end,iastep2:iastep2:end,1)))))');
axis equal;  axis tight;

[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-50:150;
intensity=zeros(length(velocities),50);
%simulate the spectrum within that box
T=1;
for t=[1:50]
    vmag=(VxC(:,:,t)'.^2+VyC(:,:,t)'.^2).^0.5;
    vmag=-1.*vmag.*sign(VxC(:,:,t)');
   pos=find([BW(1:end-1,1:end-1)]);
    vmag_small=vmag(pos);
    pos=find(isfinite(vmag_small));
    vmag_smaller=vmag_small(pos);
    v_corr(T)=mean(vmag_smaller)+std(vmag_smaller);
    T=T+1;
end

%spectra
figure(200);
hold on; plot(times, v_corr(indices),'r');
ylim([-0.2 1.3]);
%xlim([0 0.96]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;

%overlay spectra
vars={'average', 'ourVpix', 'unwrap', 'time', 'scaledMaxV', 'scaledMinV'};
load('../data/flow_phantom/b/Doppler/b_Doppler.mat', vars{:});
scaledMaxV=[scaledMaxV(end-10:end) scaledMaxV(1:end-11)];
scaledMinV=[scaledMinV(end-10:end) scaledMinV(1:end-11)];
timesDS=time;
hold on; plot(timesDS, scaledMaxV./100,'Color',[0.55 0.55 0.55]);
hold on; plot(timesDS, scaledMinV./100,'Color',[0.55 0.55 0.55]);
set(gca,'Position',[0.15 0.18 0.8 0.78]);
strings={'no corrections','with steady velocity correction','with total correction','Doppler'};
legend(strings,'Position',[0.6 0.78 0.2 0.14],'FontSize',6);
legend('boxoff')
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
set(gca,'FontSize',8);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [8.4 5]);
set(gcf, 'PaperPosition', [0 0 8.4 5]);
figfilename=['phantom_b_Delta' int2str(delta) '_corrections_Doppler_comp_2.3.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%% ----then what are the errors
v_total_error_abs=(-v_nocorr+v_corr);
v_velocity_error_abs=(-v_nocorr+v_vcorr);
v_acceleration_error_abs=(-v_vcorr+v_corr);
acc=(diff([v_corr v_corr(1)])+diff([v_corr(end) v_corr]))./(2*timestep);
error_est=v_corr./(3+v_corr).*((acc*centre(1)./(100*v_corr))-v_corr);

v_total_error=(v_nocorr-v_corr)./v_corr;
v_velocity_error=(v_nocorr-v_vcorr)./v_corr;
v_acceleration_error=(v_vcorr-v_corr)./v_corr;

figure(201);
hold on; plot(times, v_velocity_error_abs(indices),'k--');
hold on; plot(times, v_acceleration_error_abs(indices),'k:');
hold on; plot(times, v_total_error_abs(indices),'k');
% hold on; plot(times, error_est(indices),'b');
ylim([-0.3 0.3]);
strings={'steady velocity correction', 'additional acceleration correction','total correction'};
legend(strings,'Position',[0.6 0.78 0.2 0.14],'FontSize',6);
legend('boxoff')
set(gca,'FontSize',8,'Position',[0.15 0.18 0.8 0.78]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [8.4 5]);
set(gcf, 'PaperPosition', [0 0 8.4 5]);
figfilename=['phantom_b_Delta' int2str(delta) '_errors_2.3.png'];
print( '-painters', '-dpng', '-r600', figfilename);

figure(202);
hold on; plot(times, v_velocity_error(indices).*100,'k--');
hold on; plot(times, v_acceleration_error(indices).*100,'k:');
hold on; plot(times, v_total_error(indices).*100,'k');