%This script generates the images for figure 7 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/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,upperEdge,lowerEdge]=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

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=2.75e-6;
vel= womer_Vmean(timeT(pos), 20, vmean2, Rmean,nu);

flow_f(pos)'./(pi.*Rmean.^2);
vmean2=vmean./(60*1000*1000);
vel3= womer_Vmean(times, 20, vmean2, Rmean,'k',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_e_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));

%% --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;
heightBox=0.2;
widthBox=3;
[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-70:200;
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);
    for n=1:length(vmag_smaller)
        v=round(vmag_smaller(n)*100);
        pos=find(velocities==v);
        intensity(pos,T)=intensity(pos,T)+1;
    end
    averageUIV(t)=mean(vmag_smaller);
    T=T+1;
end
intensity=flipud(intensity);

%overlay spectra
vars={'average', 'ourVpix', 'unwrap', 'time', 'scaledMaxV', 'scaledMinV'};
load('../data/flow_phantom/e/Doppler/e_Doppler.mat', vars{:});
[v_long,t_long]=size(average);
velocitiesDS=(-[0:v_long]-(-64+unwrap-128)).*ourVpix.*0.01;
timesDS=time;

%spectra
fig=figure;
imagesc(timesDS, velocitiesDS, average); colormap gray;
intensity2= ind2rgb(flipud(intensity),jet);
hold on; g=imagesc(times, velocities./100,intensity2(:,indices,:));
set(g,'AlphaData',flipud(intensity(:,indices)).^0.5);
set(g,'AlphaDataMapping','scaled');
set(gca,'YDir','normal');
set(gca, 'position', [0.22 0.2 0.77 0.65],'FontSize',8,'FontWeight','normal');
hold off;
ylim([-0.7 2]);
xlim([0 1]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;
title({'Doppler spectrum      ',['overlayed with UIV spectrum \Deltat=' int2str(delta) '\tau            ']},'FontSize',8,'FontWeight','normal');

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [5.5 5]);
set(gcf, 'PaperPosition', [0 0 5.5 5]);
figfilename=['phantom_e_Delta' int2str(delta) '_Doppler_comp.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%error
background=mean(average(1,:));
sumAverageDS=(velocitiesDS(1:end-1)'*ones(1,length(timesDS))).*(average-background);
weightedAverageDS=sum(sumAverageDS)./sum(average-background);
weightedAverageDS2=interp1(timesDS,weightedAverageDS,times);
meanDiffaveDS=mean(abs(averageUIV(indices(1:end-1))-weightedAverageDS2(1:end-1)))
percentDiffaveDS=meanDiffaveDS/max(averageUIV).*100;

peak=max(average);
[a,b]=size(average);
for t=1:b
    threshDS=find(average(:,t)>(peak(t)/1.2));
    maxDS(t)=velocitiesDS(threshDS(1));
end
maxDS2=interp1(timesDS,maxDS,times);
meanDiffmaxDS=mean(abs(averageUIV(indices(1:end-1))-maxDS2(1:end-1)))
percentDiffmaxDS=meanDiffmaxDS/max(averageUIV).*100


%% -----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_e_Delta' int2str(delta) 'errors.txt'];

fid=fopen(filename,'w');
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffaveDS,percentDiffaveDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffmaxDS,percentDiffmaxDS);
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/e/UIV/2cm/conventional/phantom_e_128lines_conventional_20000depth7000000freq_out']
load([name '.mat']);

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--
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, 465, 530, 795, 855);
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);
xlim([0 1]);
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','conventional','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_e_flow_rate_comparison.png'];
print( '-painters', '-dpng', '-r600', figfilename);

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);
vel= womer_Vmean(timeT(pos), 20, vmean2, Rmean,'k',nu);

flow_f(pos)'./(pi.*Rmean.^2);
vmean2=vmean./(60*1000*1000);
vel3= womer_Vmean(times, 20, vmean2, Rmean,'k',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_e_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));

%% --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;
heightBox=0.2;
widthBox=3;
[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-70:200;
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);
    for n=1:length(vmag_smaller)
        v=round(vmag_smaller(n)*100);
        pos=find(velocities==v);
        intensity(pos,T)=intensity(pos,T)+1;
    end
    averageUIV(t)=mean(vmag_smaller);
    T=T+1;
end
intensity=flipud(intensity);

%overlay spectra
vars={'average', 'ourVpix', 'unwrap', 'time', 'scaledMaxV', 'scaledMinV'};
load('../data/flow_phantom/e/Doppler/e_Doppler.mat', vars{:});
[v_long,t_long]=size(average);
velocitiesDS=(-[0:v_long]-(-64+unwrap-128)).*ourVpix.*0.01;
timesDS=time;

%spectra
fig=figure;
imagesc(timesDS, velocitiesDS, average); colormap gray;
intensity2= ind2rgb(flipud(intensity),jet);
hold on; g=imagesc(times, velocities./100,intensity2(:,indices,:));
set(g,'AlphaData',flipud(intensity(:,indices)).^0.5);
set(g,'AlphaDataMapping','scaled');
set(gca,'YDir','normal');
set(gca, 'position', [0.22 0.2 0.77 0.65],'FontSize',8,'FontWeight','normal');
hold off;
ylim([-0.7 2]);
xlim([0 1]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;
title({'Doppler spectrum   ','overlayed with UIV spectrum conventional               '},'FontSize',8,'FontWeight','normal');

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [5.5 5]);
set(gcf, 'PaperPosition', [0 0 5.5 5]);
figfilename=['phantom_e_conventional_Doppler_comp.png'];
print( '-painters', '-dpng', '-r600', figfilename);
%error
background=mean(average(1,:));
sumAverageDS=(velocitiesDS(1:end-1)'*ones(1,length(timesDS))).*(average-background);
weightedAverageDS=sum(sumAverageDS)./sum(average-background);
weightedAverageDS2=interp1(timesDS,weightedAverageDS,times);
meanDiffaveDS=mean(abs(averageUIV(indices(1:end-1))-weightedAverageDS2(1:end-1)))
percentDiffaveDS=meanDiffaveDS/max(averageUIV).*100;

% peak=background*5;
peak=max(average);
[a,b]=size(average);
for t=1:b
    threshDS=find(average(:,t)>(peak(t)/1.2));
    maxDS(t)=velocitiesDS(threshDS(1));
end
maxDS2=interp1(timesDS,maxDS,times);
meanDiffmaxDS=mean(abs(averageUIV(indices(1:end-1))-maxDS2(1:end-1)))
percentDiffmaxDS=meanDiffmaxDS/max(averageUIV).*100

%% ------write out errors file
filename=['phantom_e_conventional_errors.txt'];
fid=fopen(filename,'w');
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffaveDS,percentDiffaveDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffmaxDS,percentDiffmaxDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanFlowerror,percentFlowerror);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanWomerror,percenterror);
fclose('all');